Commit 001f8176 authored by Praetorius, Simon's avatar Praetorius, Simon

navier stokes demo. Note: not yet working

parent 3c050d06
Pipeline #965 passed with stage
in 8 minutes and 25 seconds
......@@ -43,7 +43,6 @@ namespace AMDiS
void DirichletBC<WorldVector, Range>::finishImpl(
Matrix& matrix, VectorX& solution, VectorB& rhs,
RowBasis const& rowBasis, ColBasis const& colBasis, ValueCat)
// Dune::TypeTree::LeafNodeTag, Dune::TypeTree::LeafNodeTag)
{
using Dune::Functions::interpolate;
......
......@@ -50,7 +50,7 @@ namespace AMDiS
private: // typedefs and static constants
using GlobalBasis = typename Traits::GlobalBasis;
using GridView = typename GlobalBasis::GridView;
using Vector = DOFVectorConstView<GlobalBasis,TreePath>;
using Vector = DOFVectorConstView<Traits,TreePath>;
using Range = typename Vector::Range;
/// Dimension of the mesh
......@@ -112,10 +112,10 @@ namespace AMDiS
};
template <class Traits, class GlobalBasis, class TreePath>
template <class Traits, class TreePath>
std::shared_ptr<FileWriter<Traits,TreePath>> makeFileWriterPtr(
std::string baseName,
DOFVectorConstView<GlobalBasis,TreePath> const& dofvector)
DOFVectorConstView<Traits,TreePath> const& dofvector)
{
return std::make_shared<FileWriter<Traits,TreePath>>(baseName, dofvector);
}
......
......@@ -56,6 +56,22 @@ namespace AMDiS
, order_(order)
{}
// Copy constructor
GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
: gridFct_(that.gridFct_)
, localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
, quadCreator_(that.quadCreator_)
, order_(that.order_)
{}
// Move constructor
GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
: gridFct_(std::move(that.gridFct_))
, localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
, quadCreator_(std::move(that.quadCreator_))
, order_(std::move(that.order_))
{}
/// \brief Binds operator to `element` and `geometry`.
/**
* Binding an operator to the currently visited element in grid traversal.
......
......@@ -21,7 +21,7 @@ namespace AMDiS
// specialization for intersections
template <class I>
struct ContextImpl<I, Void_t<typename I::Entity>>
struct ContextImpl<I, Void_t<decltype(std::declval<I>().inside())>>
{
using Entity = typename I::Entity;
using Geometry = typename I::LocalGeometry;
......
......@@ -63,7 +63,7 @@ namespace AMDiS
using WorldVector = typename Element::Geometry::GlobalCoordinate;
using SystemMatrix = DOFMatrix<GlobalBasis, GlobalBasis, double>;
using SystemVector = DOFVector<GlobalBasis, double>;
using SystemVector = DOFVector<Traits, double>;
using LinearSolverType = LinearSolverInterface<typename SystemMatrix::BaseMatrix, typename SystemVector::BaseVector>;
......
......@@ -122,7 +122,7 @@ void ProblemStat<Traits>::createFileWriter()
if (!Parameters::get<std::string>(componentName + "->filename"))
return;
auto writer = makeFileWriterPtr<Traits>(componentName, this->getSolution(treePath));
auto writer = makeFileWriterPtr(componentName, this->getSolution(treePath));
filewriter.push_back(std::move(writer));
});
}
......@@ -146,7 +146,7 @@ void ProblemStat<Traits>::addMatrixOperator(
auto j = child(globalBasis->localView().tree(), makeTreePath(col));
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(op, i, j);
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i, j);
matrixOperators[i][j].element.push_back({localAssembler, factor, estFactor});
matrixOperators[i][j].changing = true;
......@@ -171,7 +171,7 @@ void ProblemStat<Traits>::addMatrixOperator(
using Intersection = typename GridView::Intersection;
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i, j);
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i, j);
matrixOperators[i][j].boundary.push_back({localAssembler, factor, estFactor, b});
matrixOperators[i][j].changing = true;
......@@ -191,7 +191,7 @@ void ProblemStat<Traits>::addVectorOperator(
auto i = child(globalBasis->localView().tree(), makeTreePath(path));
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto localAssembler = makeLocalAssemblerPtr<Element>(op, i);
auto localAssembler = makeLocalAssemblerPtr<Element>(std::move(op), i);
rhsOperators[i].element.push_back({localAssembler, factor, estFactor});
rhsOperators[i].changing = true;
......@@ -213,7 +213,7 @@ void ProblemStat<Traits>::addVectorOperator(
using Intersection = typename GridView::Intersection;
auto op = makeGridOperator(preOp, globalBasis->gridView());
auto localAssembler = makeLocalAssemblerPtr<Intersection>(op, i);
auto localAssembler = makeLocalAssemblerPtr<Intersection>(std::move(op), i);
rhsOperators[i].boundary.push_back({localAssembler, factor, estFactor, b});
rhsOperators[i].changing = true;
......
......@@ -39,24 +39,22 @@ namespace AMDiS
return v /= factor;
}
// some arithmetic operations with FieldMatrix
template <class T, int N, int M>
FieldMatrix<T,N,M> operator+(FieldMatrix<T,N,M> A, FieldMatrix<T,N,M> const& B)
template <class T>
FieldVector<T,1> operator*(FieldVector<T,1> const& v, FieldVector<T,1> const& w)
{
return A += B;
return {v[0] * w[0]};
}
template <class T, int N, int M>
FieldMatrix<T,N,M> operator-(FieldMatrix<T,N,M> A, FieldMatrix<T,N,M> const& B)
template <class T, int N>
FieldVector<T,N> operator*(FieldVector<T,1> const& factor, FieldVector<T,N> v)
{
return A -= B;
return v *= factor[0];
}
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec)
template <class T, int N>
FieldVector<T,N> operator*(FieldVector<T,N> v, FieldVector<T,1> const& factor)
{
return Dune::FMatrixHelp::mult(mat, vec);
return v *= factor[0];
}
// ----------------------------------------------------------------------------
......@@ -84,6 +82,14 @@ namespace AMDiS
return vec1.dot(vec2);
}
template <class T, int N, int M,
REQUIRES( N!=1 && M!=1 )>
auto operator*(FieldVector<T,N> v, FieldVector<T,M> const& w)
{
static_assert(M == N, "Requires vectors of the same type!");
return v.dot(w);
}
// ----------------------------------------------------------------------------
namespace Impl
......@@ -350,6 +356,9 @@ namespace AMDiS
return sqrt(dot(DF[0], DF[0]));
}
// ----------------------------------------------------------------------------
// some arithmetic operations with FieldMatrix
template <class T, int M, int N>
FieldMatrix<T,N,M> trans(FieldMatrix<T, M, N> const& A)
{
......@@ -362,6 +371,45 @@ namespace AMDiS
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator*(T scalar, FieldMatrix<T, M, N> A)
{
return A *= scalar;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator*(FieldMatrix<T, M, N> A, T scalar)
{
return A *= scalar;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator/(FieldMatrix<T, M, N> A, T scalar)
{
return A /= scalar;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator+(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A += B;
}
template <class T, int M, int N>
FieldMatrix<T,M,N> operator-(FieldMatrix<T, M, N> A, FieldMatrix<T, M, N> const& B)
{
return A -= B;
}
template <class T, int N, int M>
FieldVector<T,N> operator*(FieldMatrix<T,N,M> const& mat, FieldVector<T,M> const& vec)
{
return Dune::FMatrixHelp::mult(mat, vec);
}
template <class T, int M, int N, int L>
FieldMatrix<T,M,N> multiplies(FieldMatrix<T, M, L> const& A, FieldMatrix<T, L, N> const& B)
{
......
......@@ -80,4 +80,17 @@ namespace AMDiS
} // end namespace Category
template <class V>
constexpr bool isVector(V const&)
{
static_assert(Category::Vector<V>,"");
return Category::Vector<V>;
}
template <class V>
constexpr bool isNotVector(V const&)
{
static_assert(!Category::Vector<V>,"");
return !Category::Vector<V>;
}
} // end namespace AMDiS
......@@ -20,13 +20,13 @@ namespace AMDiS
* @{
**/
template <class GlobalBasisType, class TreePathType>
template <class Traits, class TreePathType>
class DOFVectorConstView
{
public:
using GlobalBasis = GlobalBasisType;
using GlobalBasis = typename Traits::GlobalBasis;
using TreePath = TreePathType;
using Vector = DOFVector<GlobalBasis>;
using Vector = DOFVector<Traits>;
using Tree = typename GlobalBasis::LocalView::Tree;
using SubTree = typename Dune::TypeTree::ChildForTreePath<Tree, TreePath>;
......@@ -38,7 +38,8 @@ namespace AMDiS
using Domain = typename EntitySet::GlobalCoordinate;
using Range = RangeType_t<SubTree>;
using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Range(Domain)>;
using RawSignature = typename Dune::Functions::SignatureTraits<Range(Domain)>::RawSignature;
using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
using DerivativeRange = typename DerivativeTraits::Range;
using LocalDomain = typename EntitySet::LocalCoordinate;
......@@ -177,6 +178,7 @@ namespace AMDiS
/// \brief Create a LocalFunction representing the gradient. \relates GradientLocalFunction
friend GradientLocalFunction derivative(LocalFunction const& localFunction)
{
static_assert(isValidRange<DerivativeTraits>(),"Derivative of DOFVector not defined.");
return GradientLocalFunction{*localFunction.globalFunction_};
}
......@@ -207,7 +209,7 @@ namespace AMDiS
public:
/// Constructor. Stores a pointer to the dofVector and a copy of the treePath.
DOFVectorConstView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath)
DOFVectorConstView(DOFVector<Traits> const& dofVector, TreePath const& treePath)
: dofVector_(&dofVector)
, treePath_(treePath)
, entitySet_(dofVector.getFeSpace().gridView())
......@@ -246,13 +248,13 @@ namespace AMDiS
}
/// Return const coefficient vector
DOFVector<GlobalBasis> const& coefficients() const
DOFVector<Traits> const& coefficients() const
{
return *dofVector_;
}
protected:
DOFVector<GlobalBasis> const* dofVector_;
DOFVector<Traits> const* dofVector_;
TreePath const treePath_;
EntitySet entitySet_;
......@@ -261,18 +263,18 @@ namespace AMDiS
// A mutable version of DOFVectorView
template <class GlobalBasisType, class TreePathType>
template <class Traits, class TreePathType>
class DOFVectorMutableView
: public DOFVectorConstView<GlobalBasisType, TreePathType>
: public DOFVectorConstView<Traits, TreePathType>
{
using Super = DOFVectorConstView<GlobalBasisType, TreePathType>;
using Super = DOFVectorConstView<Traits, TreePathType>;
using GlobalBasis = GlobalBasisType;
using GlobalBasis = typename Traits::GlobalBasis;
using TreePath = TreePathType;
public:
/// Constructor. Stores a pointer to the mutable `dofvector`.
DOFVectorMutableView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
DOFVectorMutableView(DOFVector<Traits>& dofVector, TreePath const& treePath)
: Super(dofVector, treePath)
, mutableDofVector_(&dofVector)
{}
......@@ -287,7 +289,7 @@ namespace AMDiS
auto&& gridFct = makeGridFunction(std::forward<Expr>(expr), basis.gridView());
DOFVector<GlobalBasis> tmp(basis, "tmp");
DOFVector<Traits> tmp(basis, "tmp");
Dune::Functions::interpolate(basis, treePath, tmp, std::forward<decltype(gridFct)>(gridFct));
// move data from temporary vector into stored DOFVector
......@@ -296,13 +298,13 @@ namespace AMDiS
}
/// Return the mutable DOFVector
DOFVector<GlobalBasis>& coefficients() { return *mutableDofVector_; }
DOFVector<Traits>& coefficients() { return *mutableDofVector_; }
/// Return the const DOFVector
using Super::coefficients;
protected:
DOFVector<GlobalBasis>* mutableDofVector_;
DOFVector<Traits>* mutableDofVector_;
};
/** @} **/
......@@ -310,34 +312,34 @@ namespace AMDiS
#ifndef DOXYGEN
// A Generator for a const \ref DOFVectorView.
template <class GlobalBasis, class TreePath>
auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector, TreePath const& treePath)
template <class Traits, class TreePath>
auto makeDOFVectorView(DOFVector<Traits> const& dofVector, TreePath const& treePath)
{
return DOFVectorConstView<GlobalBasis, TreePath>{dofVector, treePath};
return DOFVectorConstView<Traits, TreePath>{dofVector, treePath};
}
// A Generator for a mutable \ref DOFVectorView.
template <class GlobalBasis, class TreePath>
auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector, TreePath const& treePath)
template <class Traits, class TreePath>
auto makeDOFVectorView(DOFVector<Traits>& dofVector, TreePath const& treePath)
{
return DOFVectorMutableView<GlobalBasis, TreePath>{dofVector, treePath};
return DOFVectorMutableView<Traits, TreePath>{dofVector, treePath};
}
// A Generator for a const \ref DOFVectorView.
template <class GlobalBasis>
auto makeDOFVectorView(DOFVector<GlobalBasis> const& dofVector)
template <class Traits>
auto makeDOFVectorView(DOFVector<Traits> const& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DOFVectorConstView<GlobalBasis, decltype(treePath)>{dofVector, treePath};
return DOFVectorConstView<Traits, decltype(treePath)>{dofVector, treePath};
}
// A Generator for a mutable \ref DOFVectorView.
template <class GlobalBasis>
auto makeDOFVectorView(DOFVector<GlobalBasis>& dofVector)
template <class Traits>
auto makeDOFVectorView(DOFVector<Traits>& dofVector)
{
auto treePath = Dune::TypeTree::hybridTreePath();
return DOFVectorMutableView<GlobalBasis, decltype(treePath)>{dofVector, treePath};
return DOFVectorMutableView<Traits, decltype(treePath)>{dofVector, treePath};
}
#endif
......
......@@ -4,8 +4,8 @@
namespace AMDiS {
template <class GlobalBasis, class TreePath>
typename DOFVectorConstView<GlobalBasis, TreePath>::Range DOFVectorConstView<GlobalBasis, TreePath>::
template <class Traits, class TreePath>
typename DOFVectorConstView<Traits, TreePath>::Range DOFVectorConstView<Traits, TreePath>::
LocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
......@@ -58,8 +58,8 @@ LocalFunction::operator()(LocalDomain const& x) const
}
template <class GlobalBasis, class TreePath>
typename DOFVectorConstView<GlobalBasis, TreePath>::DerivativeRange DOFVectorConstView<GlobalBasis, TreePath>::
template <class Traits, class TreePath>
typename DOFVectorConstView<Traits, TreePath>::DerivativeRange DOFVectorConstView<Traits, TreePath>::
GradientLocalFunction::operator()(LocalDomain const& x) const
{
assert( bound_ );
......
......@@ -29,13 +29,16 @@ namespace AMDiS
using GridFctRange = typename GridFunction::Range;
using GridFctDomain = typename GridFunction::Domain;
using RawSignature = typename Dune::Functions::SignatureTraits<GridFctRange(GridFctDomain)>::RawSignature;
using DerivativeTraits = Dune::Functions::DefaultDerivativeTraits<RawSignature>;
using LocalFunction = std::decay_t<decltype(derivative(localFunction(std::declval<GridFunction>())))>;
public:
/// The Range of the derivative of the GridFunction
using Range = typename Dune::Functions::DefaultDerivativeTraits<RawSignature>::Range;
using Range = typename DerivativeTraits::Range;
/// The domain of the GridFunction
using Domain = GridFctDomain;
/// The EntitySet of the GridFunction
using EntitySet = typename GridFunction::EntitySet;
......@@ -43,7 +46,9 @@ namespace AMDiS
/// Constructor. Stores a copy of gridFct.
explicit DerivativeGridFunction(GridFunction const& gridFct)
: gridFct_(gridFct)
{}
{
static_assert(isValidRange<DerivativeTraits>(), "Derivative of GridFunction not defined");
}
/// NOTE: no global derivative available
Range operator()(Domain const& x) const
......
......@@ -19,16 +19,16 @@ namespace AMDiS
struct DomainType
{
using type = typename T0::Domain;
static_assert( all_of_v< std::is_same<type, typename DomainType<Ts>::type>::value... >,
"All GridFunctions must have the same Domain." );
// static_assert( all_of_v< std::is_same<type, typename DomainType<Ts>::type>::value... >,
// "All GridFunctions must have the same Domain." );
};
template <class T0, class... Ts>
struct EntitySetType
{
using type = typename T0::EntitySet;
static_assert( all_of_v< std::is_same<type, typename EntitySetType<Ts>::type>::value... >,
"All GridFunctions must have the same EntitySet." );
// static_assert( all_of_v< std::is_same<type, typename EntitySetType<Ts>::type>::value... >,
// "All GridFunctions must have the same EntitySet." );
};
} // end namespace Impl
......@@ -239,7 +239,7 @@ namespace AMDiS
{
static_assert(all_of_v<Concepts::GridFunction<GridFcts>...>,
"All passed parameters must be GridFunctions.");
static_assert(Concepts::Evaluable<Functor, typename GridFcts::Range...>,
static_assert(Concepts::Callable<Functor, typename GridFcts::Range...>,
"Range types of grid functions are not compatible with the functor.");
return FunctorGridFunction<Functor, GridFcts...>{f, gridFcts...};
}
......@@ -277,9 +277,9 @@ namespace AMDiS
* - `invokeAtQP(Operation::Plus{}, X(0), X(1));`
**/
template <class Functor, class... GridFcts>
auto invokeAtQP(Functor const& f, GridFcts const&... gridFcts)
auto invokeAtQP(Functor const& f, GridFcts&&... gridFcts)
{
return FunctorPreGridFunction<Functor, GridFcts...>{f, gridFcts...};
return FunctorPreGridFunction<Functor, std::decay_t<GridFcts>...>{f, std::forward<GridFcts>(gridFcts)...};
}
/** @} **/
......
......@@ -128,6 +128,13 @@ namespace AMDiS
} // end namespace Concepts
template <class Traits,
bool valid = !std::is_same<typename Traits::Range, Dune::Functions::InvalidRange>::value>
constexpr bool isValidRange()
{
static_assert(valid, "Traits::Range is invalid");
return valid;
}
namespace Impl
{
......
......@@ -186,7 +186,7 @@ namespace AMDiS
REQUIRES(Concepts::AnyGridFunction<Mat>)>
auto trans(Mat&& mat)
{
return invokeAtQP([](auto const& m) { return trans(m); }, std::forward<Mat>(mat));
return invokeAtQP(Operation::Trans{}, std::forward<Mat>(mat));
}
// @}
......
......@@ -17,20 +17,20 @@
namespace AMDiS
{
/// The basic container that stores a base vector and a corresponding feSpace
template <class FeSpaceType, class ValueType = double>
template <class Traits, class ValueType = double>
class DOFVector
{
using Self = DOFVector;
public:
/// The type of the \ref feSpace
using FeSpace = FeSpaceType;
using GlobalBasis = typename Traits::GlobalBasis;
/// The type of the base vector
using BaseVector = MTLDenseVector<ValueType>;
/// The index/size - type
using size_type = typename FeSpace::size_type;
using size_type = typename GlobalBasis::size_type;
/// The type of the elements of the DOFVector
using value_type = ValueType;
......@@ -39,7 +39,7 @@ namespace AMDiS
using field_type = typename BaseVector::value_type;
/// Constructor. Constructs new BaseVector.
DOFVector(FeSpace const& feSpace, std::string name)
DOFVector(GlobalBasis const& feSpace, std::string name)
: feSpace(feSpace)
, name(name)
, vector(ClonablePtr<BaseVector>::make())
......@@ -49,7 +49,7 @@ namespace AMDiS
}
/// Constructor. Takes reference to existing BaseVector
DOFVector(FeSpace const& feSpace, std::string name,
DOFVector(GlobalBasis const& feSpace, std::string name,
BaseVector& vector_ref)
: feSpace(feSpace)
, name(name)
......@@ -57,7 +57,7 @@ namespace AMDiS
{}
/// Return the basis \ref feSpace of the vector
FeSpace const& getFeSpace() const
GlobalBasis const& getFeSpace() const
{
return feSpace;
}
......@@ -154,7 +154,7 @@ namespace AMDiS
private:
/// The finite element space / basis associated with the data vector
FeSpace const& feSpace;
GlobalBasis const& feSpace;
/// The name of the DOFVector (used in file-writing)
std::string name;
......@@ -168,11 +168,11 @@ namespace AMDiS
/// Constructor a dofvector from given feSpace and name
template <class FeSpaceType, class ValueType = double>
DOFVector<FeSpaceType, ValueType>
makeDOFVector(FeSpaceType const& feSpace, std::string name)
template <class Traits, class ValueType = double>
DOFVector<Traits, ValueType>
makeDOFVector(typename Traits::GlobalBasis const& basis, std::string name)
{
return {feSpace, name};
return {basis, name};
}
} // end namespace AMDiS
......@@ -78,6 +78,23 @@ namespace AMDiS
}
};
// -------------------------------------------------------------------------
struct Trans
{
template <class M>
constexpr auto operator()(M const& mat) const
{
return trans(mat);
}
friend constexpr int order(Trans, int d)
{
return d;
}
};
/** @} **/
} // end namespace Operation
......
......@@ -38,7 +38,18 @@ namespace AMDiS
struct RangeTypeImpl<Node, false, true, false>
{
using ChildNode = typename Node::template Child<0>::type;
using type = Dune::FieldVector<RangeType_t<ChildNode>, int(Node::CHILDREN)>;
template <bool childIsLeaf, class ChildRangeType>
struct FlatType {
using type = Dune::FieldVector<ChildRangeType, int(Node::CHILDREN)>;
};
template <class T>
struct FlatType<true, Dune::FieldVector<T,1>> {
using type = Dune::FieldVector<T, int(Node::CHILDREN)>;
};
using type = typename FlatType<ChildNode::isLeaf, RangeType_t<ChildNode>>::type;
};
// Composite node
......
dimension of world: 2
stokesMesh->macro file name: ./macro/macro.stand.2d
stokesMesh->global refinements: 6
stokes->mesh: stokesMesh
stokes->solver->name: bicgstab_ell
stokes->solver->ell: 5
stokes->solver->max iteration: 1000
stokes->solver->reduction: 1e-8
stokes->solver->restart: 50
stokes->solver->info: 2
stokes->output[0]->filename: stokes_u.2d
stokes->output[0]->name: u
stokes->output[0]->subsampling: 0
stokes->output[0]->output directory: ./output
stokes->output[0]->ParaView format: 1
stokes->output[0]->ParaView mode: 1
stokes->output[1]->filename: stokes_p.2d
stokes->output[1]->name: p
stokes->output[1]->output directory: ./output
stokes->output[1]->ParaView format: 1
stokes->output[1]->ParaView mode: 1
stokes->viscosity: 0.1
stokes->viscosity: 1.0
stokes->boundary velocity: 10.0
adapt->timestep: 0.1
adapt->start time: 0.0
adapt->end time: 10.0
......@@ -6,7 +6,8 @@
#include <dune/amdis/AdaptInstationary.hpp>
#include <dune/amdis/ProblemStat.hpp>
#include <dune/amdis/ProblemInstat.hpp>
#include <dune/amdis/Terms.hpp>
#include <dune/amdis/Operators.hpp>
#include <dune/amdis/assembler/StokesOperator.hpp>
#ifdef DOW
#undef DOW
......@@ -19,12 +20,14 @@
using namespace AMDiS;
using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
using StokesParam = TaylorHoodBasis<Grid::LeafGridView>;
using StokesProblem = ProblemStat<StokesParam>;
using StokesProblemInstat = ProblemInstat<StokesParam>;
struct NavierStokesBasis
{
using Grid = Dune::AlbertaGrid<AMDIS_DIM,AMDIS_DOW>;
using GlobalBasis = typename TaylorHoodBasis<Grid::LeafGridView>::GlobalBasis;
};
using Op = StokesProblem::ElementOperator;
using StokesProblem = ProblemStat<NavierStokesBasis>;
using StokesProblemInstat = ProblemInstat<NavierStokesBasis>;
int main(int argc, char** argv)
{
......@@ -57,20 +60,10 @@ int main(int argc, char** argv)
auto opTimeOld = makeOperator(tag::testvec{}, density * prob.getSolution(_v));
prob.addVectorOperator(opTimeOld, _v);
// sum_i <grad(u_i),grad(v_i)>
for (std::size_t i = 0; i < DOW; ++i) {
auto opL = makeOperator(tag::gradtest_gradtrial{}, viscosity);
prob.addMatrixOperator(opL, treepath(_v,i), treepath(_v,i));