Skip to content
Snippets Groups Projects
Commit 1e8b463d authored by Elias Pipping's avatar Elias Pipping
Browse files

[Cleanup] Drop Type suffix where possible

parent 3c0fcc10
Branches
No related tags found
No related merge requests found
Showing
with 298 additions and 326 deletions
......@@ -16,10 +16,10 @@ template <size_t dim> class EllipticEnergy {
using SmallVector = FieldVector<double, dim>;
using SmallMatrix = FieldMatrix<double, dim, dim>;
using NonlinearityType = LocalFriction<dim>;
using Nonlinearity = LocalFriction<dim>;
EllipticEnergy(SmallMatrix const &A, SmallVector const &b,
shared_ptr<NonlinearityType const> phi, size_t ignore = dim)
shared_ptr<Nonlinearity const> phi, size_t ignore = dim)
: A(A), b(b), phi(phi), ignore(ignore) {}
double operator()(SmallVector const &v) const {
......@@ -28,7 +28,7 @@ template <size_t dim> class EllipticEnergy {
SmallMatrix const &A;
SmallVector const &b;
shared_ptr<NonlinearityType const> const phi;
shared_ptr<Nonlinearity const> const phi;
/* Dimension that should be ignored; goes from 0 to dim-1; the
special value dim means that no dimension should be ignored */
size_t const ignore;
......
......@@ -12,21 +12,20 @@
#include "localfriction.hh"
namespace Dune {
template <class MatrixTypeTEMPLATE, class VectorTypeTEMPLATE>
class GlobalNonlinearity {
template <class Matrix, class Vector> class GlobalNonlinearity {
protected:
using SingletonVectorType = BlockVector<FieldVector<double, 1>>;
using SingletonVector = BlockVector<FieldVector<double, 1>>;
public:
using IndexSet = Dune::MatrixIndexSet;
using MatrixType = MatrixTypeTEMPLATE;
using VectorType = VectorTypeTEMPLATE;
using LocalMatrixType = typename MatrixType::block_type;
using LocalVectorType = typename VectorType::block_type;
using MatrixType = Matrix;
using VectorType = Vector;
using LocalMatrix = typename Matrix::block_type;
using LocalVectorType = typename Vector::block_type;
size_t static const block_size = LocalVectorType::dimension;
using FrictionType = LocalFriction<block_size>;
using Friction = LocalFriction<block_size>;
double operator()(VectorType const &x) const {
double operator()(Vector const &x) const {
double tmp = 0;
for (size_t i = 0; i < x.size(); ++i) {
auto const res = restriction(i);
......@@ -40,20 +39,20 @@ class GlobalNonlinearity {
*/
virtual shared_ptr<LocalFriction<block_size>> restriction(size_t i) const = 0;
void addHessian(VectorType const &v, MatrixType &hessian) const {
void addHessian(Vector const &v, Matrix &hessian) const {
for (size_t i = 0; i < v.size(); ++i) {
auto const res = restriction(i);
res->addHessian(v[i], hessian[i][i]);
}
}
void directionalDomain(VectorType const &, VectorType const &,
void directionalDomain(Vector const &, Vector const &,
Interval<double> &dom) const {
dom[0] = -std::numeric_limits<double>::max();
dom[1] = std::numeric_limits<double>::max();
}
void directionalSubDiff(VectorType const &u, VectorType const &v,
void directionalSubDiff(Vector const &u, Vector const &v,
Interval<double> &subdifferential) const {
subdifferential[0] = subdifferential[1] = 0;
for (size_t i = 0; i < u.size(); ++i) {
......@@ -70,19 +69,19 @@ class GlobalNonlinearity {
indices.add(i, i);
}
void addGradient(VectorType const &v, VectorType &gradient) const {
void addGradient(Vector const &v, Vector &gradient) const {
for (size_t i = 0; i < v.size(); ++i) {
auto const res = restriction(i);
res->addGradient(v[i], gradient[i]);
}
}
double regularity(size_t i, typename VectorType::block_type const &x) const {
double regularity(size_t i, typename Vector::block_type const &x) const {
auto const res = restriction(i);
return res->regularity(x);
}
virtual void updateLogState(SingletonVectorType const &logState) = 0;
virtual void updateLogState(SingletonVector const &logState) = 0;
};
}
#endif
......@@ -15,34 +15,32 @@
#include "frictionpotential.hh"
namespace Dune {
template <class MatrixType, class VectorType>
class GlobalRuinaNonlinearity
: public GlobalNonlinearity<MatrixType, VectorType> {
template <class Matrix, class Vector>
class GlobalRuinaNonlinearity : public GlobalNonlinearity<Matrix, Vector> {
public:
using GlobalNonlinearity<MatrixType, VectorType>::block_size;
using typename GlobalNonlinearity<MatrixType, VectorType>::FrictionType;
using GlobalNonlinearity<Matrix, Vector>::block_size;
using typename GlobalNonlinearity<Matrix, Vector>::Friction;
private:
using typename GlobalNonlinearity<MatrixType,
VectorType>::SingletonVectorType;
using typename GlobalNonlinearity<Matrix, Vector>::SingletonVector;
public:
GlobalRuinaNonlinearity(Dune::BitSetVector<1> const &frictionalNodes,
SingletonVectorType const &nodalIntegrals,
SingletonVector const &nodalIntegrals,
FrictionData const &fd)
: restrictions(nodalIntegrals.size()) {
auto trivialNonlinearity =
make_shared<FrictionType>(make_shared<TrivialFunction>());
make_shared<Friction>(make_shared<TrivialFunction>());
for (size_t i = 0; i < restrictions.size(); ++i) {
restrictions[i] =
frictionalNodes[i][0]
? make_shared<FrictionType>(
? make_shared<Friction>(
make_shared<FrictionPotential>(nodalIntegrals[i], fd))
: trivialNonlinearity;
}
}
void updateLogState(SingletonVectorType const &logState) override {
void updateLogState(SingletonVector const &logState) override {
for (size_t i = 0; i < restrictions.size(); ++i)
restrictions[i]->updateLogState(logState[i]);
}
......@@ -50,12 +48,12 @@ class GlobalRuinaNonlinearity
/*
Return a restriction of the outer function to the i'th node.
*/
shared_ptr<FrictionType> restriction(size_t i) const override {
shared_ptr<Friction> restriction(size_t i) const override {
return restrictions[i];
}
private:
std::vector<shared_ptr<FrictionType>> restrictions;
std::vector<shared_ptr<Friction>> restrictions;
};
}
#endif
......@@ -14,8 +14,8 @@
#include "frictionpotential.hh"
// In order to compute (x * y) / |x|, compute x/|x| first
template <class VectorType>
double dotFirstNormalised(VectorType const &x, VectorType const &y) {
template <class Vector>
double dotFirstNormalised(Vector const &x, Vector const &y) {
double const xnorm = x.two_norm();
if (xnorm <= 0.0)
return 0.0;
......
......@@ -18,9 +18,9 @@ void descentMinimisation(Functional const &J,
typename Functional::SmallVector const &v,
Bisection const &bisection) {
using SmallVector = typename Functional::SmallVector;
using LocalNonlinearityType = typename Functional::NonlinearityType;
using LocalNonlinearity = typename Functional::Nonlinearity;
MyDirectionalConvexFunction<LocalNonlinearityType> const JRest(
MyDirectionalConvexFunction<LocalNonlinearity> const JRest(
computeDirectionalA(J.A, v), computeDirectionalb(J.A, J.b, x, v), *J.phi,
x, v);
// }}}
......
......@@ -19,23 +19,22 @@
#include "ellipticenergy.hh"
/* Just for debugging */
template <class MatrixType, class VectorType>
double computeEnergy(
MatrixType const &A, VectorType const &x, VectorType const &b,
Dune::GlobalNonlinearity<MatrixType, VectorType> const &phi) {
template <class Matrix, class Vector>
double computeEnergy(Matrix const &A, Vector const &x, Vector const &b,
Dune::GlobalNonlinearity<Matrix, Vector> const &phi) {
return computeEnergy(A, x, b) + phi(x);
}
/** \brief Base class for problems where each block can be solved with a
* modified gradient method */
template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
template <class ConvexProblem> class MyBlockProblem {
public:
using ConvexProblemType = ConvexProblemTypeTEMPLATE;
using VectorType = typename ConvexProblemType::VectorType;
using MatrixType = typename ConvexProblemType::MatrixType;
using LocalVectorType = typename ConvexProblemType::LocalVectorType;
using LocalMatrixType = typename ConvexProblemType::LocalMatrixType;
size_t static const block_size = ConvexProblemType::block_size;
using ConvexProblemType = ConvexProblem;
using VectorType = typename ConvexProblem::VectorType;
using MatrixType = typename ConvexProblem::MatrixType;
using LocalVector = typename ConvexProblem::LocalVectorType;
using LocalMatrix = typename ConvexProblem::LocalMatrixType;
size_t static const block_size = ConvexProblem::block_size;
size_t static const coarse_block_size = block_size;
/** \brief Solves one local system using a modified gradient method */
......@@ -44,10 +43,8 @@ template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
struct Linearization {
size_t static const block_size = coarse_block_size;
using LocalMatrixType =
typename MyBlockProblem<ConvexProblemType>::LocalMatrixType;
using MatrixType =
Dune::BCRSMatrix<typename Linearization::LocalMatrixType>;
using LocalMatrix = typename MyBlockProblem<ConvexProblem>::LocalMatrix;
using MatrixType = Dune::BCRSMatrix<typename Linearization::LocalMatrix>;
using VectorType =
Dune::BlockVector<Dune::FieldVector<double, Linearization::block_size>>;
using BitVectorType = Dune::BitSetVector<Linearization::block_size>;
......@@ -60,7 +57,7 @@ template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
};
MyBlockProblem(Dune::ParameterTree const &parset,
ConvexProblemType const &problem)
ConvexProblem const &problem)
: parset(parset), problem(problem) {
bisection = Bisection();
}
......@@ -196,7 +193,7 @@ template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
Dune::ParameterTree const &parset;
// problem data
ConvexProblemType const &problem;
ConvexProblem const &problem;
// commonly used minimization stuff
Bisection bisection;
......@@ -205,8 +202,8 @@ template <class ConvexProblemTypeTEMPLATE> class MyBlockProblem {
};
/** \brief Solves one local system using a scalar Gauss-Seidel method */
template <class ConvexProblemTypeTEMPLATE>
class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
template <class ConvexProblem>
class MyBlockProblem<ConvexProblem>::IterateObject {
friend class MyBlockProblem;
protected:
......@@ -215,7 +212,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
* \param problem The problem including quadratic part and nonlinear part
*/
IterateObject(Dune::ParameterTree const &parset, Bisection const &bisection,
ConvexProblemType const &problem)
ConvexProblem const &problem)
: problem(problem),
bisection(bisection),
localsteps(parset.get<size_t>("localsolver.steps")) {}
......@@ -228,7 +225,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
}
/** \brief Update the i-th block of the current iterate */
void updateIterate(LocalVectorType const &ui, size_t i) {
void updateIterate(LocalVector const &ui, size_t i) {
u[i] = ui;
return;
}
......@@ -239,7 +236,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
* \param ignore Set of degrees of freedom to leave untouched
*/
void solveLocalProblem(
LocalVectorType &ui, size_t m,
LocalVector &ui, size_t m,
typename Dune::BitSetVector<block_size>::const_reference ignore) {
{
size_t ic =
......@@ -258,8 +255,8 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
assert(false);
}
LocalMatrixType const *localA = nullptr;
LocalVectorType localb(problem.f[m]);
LocalMatrix const *localA = nullptr;
LocalVector localb(problem.f[m]);
auto const end = problem.A[m].end();
for (auto it = problem.A[m].begin(); it != end; ++it) {
......@@ -279,7 +276,7 @@ class MyBlockProblem<ConvexProblemTypeTEMPLATE>::IterateObject {
private:
// problem data
ConvexProblemType const &problem;
ConvexProblem const &problem;
// commonly used minimization stuff
Bisection bisection;
......
......@@ -15,24 +15,24 @@
localb = <b - Au, v>
*/
template <class MatrixType, class VectorType>
double computeDirectionalA(MatrixType const &A, VectorType const &v) {
template <class Matrix, class Vector>
double computeDirectionalA(Matrix const &A, Vector const &v) {
return Arithmetic::Axy(A, v, v);
}
template <class MatrixType, class VectorType>
double computeDirectionalb(MatrixType const &A, VectorType const &b,
VectorType const &u, VectorType const &v) {
template <class Matrix, class Vector>
double computeDirectionalb(Matrix const &A, Vector const &b, Vector const &u,
Vector const &v) {
return Arithmetic::bmAxy(A, b, u, v);
}
template <class NonlinearityType> class MyDirectionalConvexFunction {
template <class Nonlinearity> class MyDirectionalConvexFunction {
public:
using VectorType = typename NonlinearityType::VectorType;
using MatrixType = typename NonlinearityType::MatrixType;
using Vector = typename Nonlinearity::VectorType;
using Matrix = typename Nonlinearity::MatrixType;
MyDirectionalConvexFunction(double A, double b, NonlinearityType const &phi,
VectorType const &u, VectorType const &v)
MyDirectionalConvexFunction(double A, double b, Nonlinearity const &phi,
Vector const &u, Vector const &v)
: A(A), b(b), phi(phi), u(u), v(v) {
phi.directionalDomain(u, v, dom);
}
......@@ -42,7 +42,7 @@ template <class NonlinearityType> class MyDirectionalConvexFunction {
double linearPart() const { return b; }
void subDiff(double x, Interval<double> &D) const {
VectorType uxv = u;
Vector uxv = u;
Arithmetic::addProduct(uxv, x, v);
phi.directionalSubDiff(uxv, v, D);
D[0] += A * x - b;
......@@ -58,9 +58,9 @@ template <class NonlinearityType> class MyDirectionalConvexFunction {
double b;
private:
NonlinearityType const &phi;
VectorType const &u;
VectorType const &v;
Nonlinearity const &phi;
Vector const &u;
Vector const &v;
Interval<double> dom;
};
......
......@@ -11,33 +11,31 @@
#include "assemblers.hh"
// Assembles Neumann boundary term in f
template <class GridView, class LocalVectorType, class AssemblerType>
void assembleNeumann(GridView const &gridView, AssemblerType const &assembler,
template <class GridView, class LocalVector, class Assembler>
void assembleNeumann(GridView const &gridView, Assembler const &assembler,
Dune::BitSetVector<1> const &neumannNodes,
Dune::BlockVector<LocalVectorType> &f,
Dune::BlockVector<LocalVector> &f,
Dune::VirtualFunction<double, double> const &neumann,
double relativeTime) { // constant sample function on
// neumann boundary
BoundaryPatch<GridView> const neumannBoundary(gridView, neumannNodes);
LocalVectorType localNeumann(0);
LocalVector localNeumann(0);
neumann.evaluate(relativeTime, localNeumann[0]);
ConstantFunction<LocalVectorType, LocalVectorType> const fNeumann(
localNeumann);
NeumannBoundaryAssembler<typename GridView::Grid, LocalVectorType>
ConstantFunction<LocalVector, LocalVector> const fNeumann(localNeumann);
NeumannBoundaryAssembler<typename GridView::Grid, LocalVector>
neumannBoundaryAssembler(fNeumann);
assembler.assembleBoundaryFunctional(neumannBoundaryAssembler, f,
neumannBoundary);
}
// Assembles constant 1-function on frictional boundary in nodalIntegrals
template <class GridView, class LocalVectorType, class AssemblerType>
template <class GridView, class LocalVector, class Assembler>
Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
assembleFrictionWeightsal(GridView const &gridView,
AssemblerType const &assembler,
assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
Dune::BitSetVector<1> const &frictionalNodes) {
using Singleton = Dune::FieldVector<double, 1>;
BoundaryPatch<GridView> const frictionalBoundary(gridView, frictionalNodes);
ConstantFunction<LocalVectorType, Singleton> const constantOneFunction(1);
ConstantFunction<LocalVector, Singleton> const constantOneFunction(1);
NeumannBoundaryAssembler<typename GridView::Grid, Singleton>
frictionalBoundaryAssembler(constantOneFunction);
......@@ -47,14 +45,12 @@ assembleFrictionWeightsal(GridView const &gridView,
return nodalIntegrals;
}
template <class MatrixType, class VectorType>
Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assembleNonlinearity(
template <class Matrix, class Vector>
Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
FrictionData const &fd) {
return Dune::make_shared<
Dune::GlobalRuinaNonlinearity<MatrixType, VectorType>>(
return Dune::make_shared<Dune::GlobalRuinaNonlinearity<Matrix, Vector>>(
frictionalNodes, nodalIntegrals, fd);
}
......
......@@ -11,22 +11,20 @@
#include <dune/tectonic/globalnonlinearity.hh>
template <class GridView, class LocalVectorType, class AssemblerType>
void assembleNeumann(GridView const &gridView, AssemblerType const &assembler,
template <class GridView, class LocalVector, class Assembler>
void assembleNeumann(GridView const &gridView, Assembler const &assembler,
Dune::BitSetVector<1> const &neumannNodes,
Dune::BlockVector<LocalVectorType> &f,
Dune::BlockVector<LocalVector> &f,
Dune::VirtualFunction<double, double> const &neumann,
double relativeTime);
template <class GridView, class LocalVectorType, class AssemblerType>
template <class GridView, class LocalVector, class Assembler>
Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
assembleFrictionWeightsal(GridView const &gridView,
AssemblerType const &assembler,
assembleFrictionWeightsal(GridView const &gridView, Assembler const &assembler,
Dune::BitSetVector<1> const &frictionalNodes);
template <class MatrixType, class VectorType>
Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assembleNonlinearity(
template <class Matrix, class Vector>
Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>> assembleNonlinearity(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
FrictionData const &fd);
......
......@@ -15,27 +15,27 @@
using SmallVector = Dune::FieldVector<double, DIM>;
using SmallMatrix = Dune::FieldMatrix<double, DIM, DIM>;
using MatrixType = Dune::BCRSMatrix<SmallMatrix>;
using VectorType = Dune::BlockVector<SmallVector>;
using Matrix = Dune::BCRSMatrix<SmallMatrix>;
using Vector = Dune::BlockVector<SmallVector>;
using GridType = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
using GridView = GridType::LeafGridView;
using Grid = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
using GridView = Grid::LeafGridView;
using P1Basis = P1NodalBasis<GridView, double>;
using AssemblerType = Assembler<P1Basis, P1Basis>;
using MyAssembler = Assembler<P1Basis, P1Basis>;
template void assembleNeumann<GridView, SmallVector, AssemblerType>(
GridView const &gridView, AssemblerType const &assembler,
template void assembleNeumann<GridView, SmallVector, MyAssembler>(
GridView const &gridView, MyAssembler const &assembler,
Dune::BitSetVector<1> const &neumannNodes,
Dune::BlockVector<SmallVector> &f,
Dune::VirtualFunction<double, double> const &neumann, double relativeTime);
template Dune::shared_ptr<Dune::BlockVector<Dune::FieldVector<double, 1>>>
assembleFrictionWeightsal<GridView, SmallVector, AssemblerType>(
GridView const &gridView, AssemblerType const &assembler,
assembleFrictionWeightsal<GridView, SmallVector, MyAssembler>(
GridView const &gridView, MyAssembler const &assembler,
Dune::BitSetVector<1> const &frictionalNodes);
template Dune::shared_ptr<Dune::GlobalNonlinearity<MatrixType, VectorType>>
assembleNonlinearity<MatrixType, VectorType>(
template Dune::shared_ptr<Dune::GlobalNonlinearity<Matrix, Vector>>
assembleNonlinearity<Matrix, Vector>(
Dune::BitSetVector<1> const &frictionalNodes,
Dune::BlockVector<Dune::FieldVector<double, 1>> const &nodalIntegrals,
FrictionData const &fd);
......@@ -3,18 +3,18 @@
#include <dune/common/exceptions.hh>
#include <dune/common/typetraits.hh>
template <class EnumType> struct StringToEnum : public Dune::NotImplemented {};
template <class Enum> struct StringToEnum : public Dune::NotImplemented {};
template <class EnumType>
template <class Enum>
typename Dune::enable_if<
!Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<EnumType>>::value,
!Dune::IsBaseOf<Dune::NotImplemented, StringToEnum<Enum>>::value,
std::istream &>::type
operator>>(std::istream &lhs, EnumType &e) {
operator>>(std::istream &lhs, Enum &e) {
std::string s;
lhs >> s;
try {
e = StringToEnum<EnumType>::convert(s);
e = StringToEnum<Enum>::convert(s);
}
catch (typename Dune::Exception) {
lhs.setstate(std::ios_base::failbit);
......
......@@ -10,9 +10,9 @@ double computeCOF(FrictionData const &fd, double V, double logState) {
return (mu > 0) ? mu : 0;
}
template <class BitVectorType>
FrictionWriter<BitVectorType>::FrictionWriter(
FrictionData const &_fd, BitVectorType const &_boundaryNodes)
template <class BitVector>
FrictionWriter<BitVector>::FrictionWriter(FrictionData const &_fd,
BitVector const &_boundaryNodes)
: coefficientWriter("coefficients", std::fstream::out),
displacementWriter("displacements", std::fstream::out),
stateWriter("states", std::fstream::out),
......@@ -20,19 +20,17 @@ FrictionWriter<BitVectorType>::FrictionWriter(
fd(_fd),
boundaryNodes(_boundaryNodes) {}
template <class BitVectorType>
FrictionWriter<BitVectorType>::~FrictionWriter() {
template <class BitVector> FrictionWriter<BitVector>::~FrictionWriter() {
stateWriter.close();
displacementWriter.close();
velocityWriter.close();
coefficientWriter.close();
}
template <class BitVectorType>
template <class SingletonVectorType, class VectorType>
void FrictionWriter<BitVectorType>::writeInfo(SingletonVectorType const &alpha,
VectorType const &u,
VectorType const &v) {
template <class BitVector>
template <class SingletonVector, class Vector>
void FrictionWriter<BitVector>::writeInfo(SingletonVector const &alpha,
Vector const &u, Vector const &v) {
for (size_t i = 0; i < boundaryNodes.size(); ++i) {
if (!boundaryNodes[i][0])
continue;
......
......@@ -5,15 +5,15 @@
#include <dune/tectonic/frictiondata.hh>
template <class BitVectorType> class FrictionWriter {
template <class BitVector> class FrictionWriter {
public:
FrictionWriter(FrictionData const &_fd, BitVectorType const &_boundaryNodes);
FrictionWriter(FrictionData const &_fd, BitVector const &_boundaryNodes);
~FrictionWriter();
template <class SingletonVectorType, class VectorType>
void writeInfo(SingletonVectorType const &alpha, VectorType const &u,
VectorType const &v);
template <class SingletonVector, class Vector>
void writeInfo(SingletonVector const &alpha, Vector const &u,
Vector const &v);
private:
std::fstream coefficientWriter;
......@@ -21,6 +21,6 @@ template <class BitVectorType> class FrictionWriter {
std::fstream stateWriter;
std::fstream velocityWriter;
FrictionData const &fd;
BitVectorType const &boundaryNodes;
BitVector const &boundaryNodes;
};
#endif
#include <dune/common/bitsetvector.hh>
#include <dune/istl/bvector.hh>
using BitVectorType = Dune::BitSetVector<1>;
using SingletonVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using VectorType2 = Dune::BlockVector<Dune::FieldVector<double, 2>>;
using VectorType3 = Dune::BlockVector<Dune::FieldVector<double, 3>>;
using BitVector = Dune::BitSetVector<1>;
using SingletonVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using Vector2 = Dune::BlockVector<Dune::FieldVector<double, 2>>;
using Vector3 = Dune::BlockVector<Dune::FieldVector<double, 3>>;
template class FrictionWriter<BitVectorType>;
template class FrictionWriter<BitVector>;
template void FrictionWriter<BitVectorType>::writeInfo(
SingletonVectorType const &alpha, VectorType2 const &u,
VectorType2 const &v);
template void FrictionWriter<BitVectorType>::writeInfo(
SingletonVectorType const &alpha, VectorType3 const &u,
VectorType3 const &v);
template void FrictionWriter<BitVector>::writeInfo(SingletonVector const &alpha,
Vector2 const &u,
Vector2 const &v);
template void FrictionWriter<BitVector>::writeInfo(SingletonVector const &alpha,
Vector3 const &u,
Vector3 const &v);
......@@ -100,24 +100,21 @@
size_t const dims = DIM;
template <class VectorType, class MatrixType, class FunctionType, int dimension>
Dune::shared_ptr<
TimeSteppingScheme<VectorType, MatrixType, FunctionType, dimension>>
template <class Vector, class Matrix, class Function, int dimension>
Dune::shared_ptr<TimeSteppingScheme<Vector, Matrix, Function, dimension>>
initTimeStepper(Config::scheme scheme,
FunctionType const &velocityDirichletFunction,
Function const &velocityDirichletFunction,
Dune::BitSetVector<dimension> const &velocityDirichletNodes,
MatrixType const &massMatrix, MatrixType const &stiffnessMatrix,
MatrixType const &dampingMatrix, VectorType const &u_initial,
VectorType const &v_initial, VectorType const &a_initial) {
Matrix const &massMatrix, Matrix const &stiffnessMatrix,
Matrix const &dampingMatrix, Vector const &u_initial,
Vector const &v_initial, Vector const &a_initial) {
switch (scheme) {
case Config::Newmark:
return Dune::make_shared<
Newmark<VectorType, MatrixType, FunctionType, dims>>(
return Dune::make_shared<Newmark<Vector, Matrix, Function, dims>>(
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
a_initial, velocityDirichletNodes, velocityDirichletFunction);
case Config::BackwardEuler:
return Dune::make_shared<
BackwardEuler<VectorType, MatrixType, FunctionType, dims>>(
return Dune::make_shared<BackwardEuler<Vector, Matrix, Function, dims>>(
stiffnessMatrix, massMatrix, dampingMatrix, u_initial, v_initial,
velocityDirichletNodes, velocityDirichletFunction);
default:
......@@ -125,20 +122,16 @@ initTimeStepper(Config::scheme scheme,
}
}
template <class SingletonVectorType, class VectorType>
Dune::shared_ptr<StateUpdater<SingletonVectorType, VectorType>>
initStateUpdater(Config::stateModel model,
SingletonVectorType const &alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes,
FrictionData const &fd) {
template <class SingletonVector, class Vector>
Dune::shared_ptr<StateUpdater<SingletonVector, Vector>> initStateUpdater(
Config::stateModel model, SingletonVector const &alpha_initial,
Dune::BitSetVector<1> const &frictionalNodes, FrictionData const &fd) {
switch (model) {
case Config::Dieterich:
return Dune::make_shared<
DieterichStateUpdater<SingletonVectorType, VectorType>>(
return Dune::make_shared<DieterichStateUpdater<SingletonVector, Vector>>(
alpha_initial, frictionalNodes, fd.L);
case Config::Ruina:
return Dune::make_shared<
RuinaStateUpdater<SingletonVectorType, VectorType>>(
return Dune::make_shared<RuinaStateUpdater<SingletonVector, Vector>>(
alpha_initial, frictionalNodes, fd.L);
default:
assert(false);
......@@ -162,32 +155,31 @@ int main(int argc, char *argv[]) {
using SmallVector = Dune::FieldVector<double, dims>;
using SmallMatrix = Dune::FieldMatrix<double, dims, dims>;
using SmallSingletonMatrix = Dune::FieldMatrix<double, 1, 1>;
using MatrixType = Dune::BCRSMatrix<SmallMatrix>;
using SingletonMatrixType = Dune::BCRSMatrix<SmallSingletonMatrix>;
using VectorType = Dune::BlockVector<SmallVector>;
using SingletonVectorType = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using NonlinearityType = Dune::GlobalNonlinearity<MatrixType, VectorType>;
using Matrix = Dune::BCRSMatrix<SmallMatrix>;
using SingletonMatrix = Dune::BCRSMatrix<SmallSingletonMatrix>;
using Vector = Dune::BlockVector<SmallVector>;
using SingletonVector = Dune::BlockVector<Dune::FieldVector<double, 1>>;
using Nonlinearity = Dune::GlobalNonlinearity<Matrix, Vector>;
auto const E = parset.get<double>("body.E");
auto const nu = parset.get<double>("body.nu");
// {{{ Set up grid
using GridType =
Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
using Grid = Dune::ALUGrid<dims, dims, Dune::simplex, Dune::nonconforming>;
Dune::FieldVector<typename GridType::ctype, dims> lowerLeft(0);
Dune::FieldVector<typename GridType::ctype, dims> upperRight(1);
Dune::FieldVector<typename Grid::ctype, dims> lowerLeft(0);
Dune::FieldVector<typename Grid::ctype, dims> upperRight(1);
upperRight[0] = parset.get<size_t>("body.width");
upperRight[1] = parset.get<size_t>("body.height");
Dune::shared_ptr<GridType> grid;
Dune::shared_ptr<Grid> grid;
{
Dune::array<unsigned int, dims> elements;
std::fill(elements.begin(), elements.end(), 1);
elements[0] = parset.get<size_t>("body.width");
elements[1] = parset.get<size_t>("body.height");
grid = Dune::StructuredGridFactory<GridType>::createSimplexGrid(
grid = Dune::StructuredGridFactory<Grid>::createSimplexGrid(
lowerLeft, upperRight, elements);
}
......@@ -195,7 +187,7 @@ int main(int argc, char *argv[]) {
grid->globalRefine(refinements);
size_t const fineVertexCount = grid->size(grid->maxLevel(), dims);
using GridView = GridType::LeafGridView;
using GridView = Grid::LeafGridView;
GridView const leafView = grid->leafView();
// }}}
......@@ -215,7 +207,7 @@ int main(int argc, char *argv[]) {
Dune::BitSetVector<1> neumannNodes(fineVertexCount, false);
Dune::BitSetVector<1> frictionalNodes(fineVertexCount, false);
VectorType vertexCoordinates(fineVertexCount);
Vector vertexCoordinates(fineVertexCount);
{
Dune::MultipleCodimMultipleGeomTypeMapper<
GridView, Dune::MCMGVertexLayout> const vertexMapper(leafView);
......@@ -247,8 +239,8 @@ int main(int argc, char *argv[]) {
BoundaryPatch<GridView> const frictionalBoundary(leafView, frictionalNodes);
// Set up functions for time-dependent boundary conditions
using FunctionType = Dune::VirtualFunction<double, double>;
using FunctionMap = SharedPointerMap<std::string, FunctionType>;
using Function = Dune::VirtualFunction<double, double>;
using FunctionMap = SharedPointerMap<std::string, Function>;
FunctionMap functions;
{
initPython();
......@@ -262,14 +254,14 @@ int main(int argc, char *argv[]) {
// Set up normal stress, mass matrix, and gravity functional
double normalStress;
MatrixType M;
EnergyNorm<MatrixType, VectorType> const MNorm(M);
VectorType gravityFunctional;
Matrix M;
EnergyNorm<Matrix, Vector> const MNorm(M);
Vector gravityFunctional;
{
double const gravity = 9.81;
double const density = parset.get<double>("body.density");
{
MassAssembler<GridType, P1Basis::LocalFiniteElement,
MassAssembler<Grid, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement,
Dune::ScaledIdentityMatrix<double, dims>> const localMass;
p1Assembler.assembleOperator(localMass, M);
......@@ -295,8 +287,7 @@ int main(int argc, char *argv[]) {
weightedGravitationalDirection[1] = -density * gravity;
ConstantFunction<SmallVector, SmallVector> const gravityFunction(
weightedGravitationalDirection);
L2FunctionalAssembler<GridType, P1Basis::LocalFiniteElement,
SmallVector>
L2FunctionalAssembler<Grid, P1Basis::LocalFiniteElement, SmallVector>
gravityFunctionalAssembler(gravityFunction);
p1Assembler.assembleFunctional(gravityFunctionalAssembler,
gravityFunctional);
......@@ -305,80 +296,80 @@ int main(int argc, char *argv[]) {
FrictionData const frictionData(parset.sub("boundary.friction"),
normalStress);
MatrixType A;
EnergyNorm<MatrixType, VectorType> const ANorm(A);
Matrix A;
EnergyNorm<Matrix, Vector> const ANorm(A);
{
StVenantKirchhoffAssembler<GridType, P1Basis::LocalFiniteElement,
StVenantKirchhoffAssembler<Grid, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const
localStiffness(E, nu);
p1Assembler.assembleOperator(localStiffness, A);
}
MatrixType C;
Matrix C;
{
ViscosityAssembler<GridType, P1Basis::LocalFiniteElement,
ViscosityAssembler<Grid, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement> const
localViscosity(parset.get<double>("body.shearViscosity"),
parset.get<double>("body.bulkViscosity"));
p1Assembler.assembleOperator(localViscosity, C);
}
SingletonMatrixType frictionalBoundaryMassMatrix;
EnergyNorm<SingletonMatrixType, SingletonVectorType> const stateEnergyNorm(
SingletonMatrix frictionalBoundaryMassMatrix;
EnergyNorm<SingletonMatrix, SingletonVector> const stateEnergyNorm(
frictionalBoundaryMassMatrix);
{
BoundaryMassAssembler<
GridType, BoundaryPatch<GridView>, P1Basis::LocalFiniteElement,
Grid, BoundaryPatch<GridView>, P1Basis::LocalFiniteElement,
P1Basis::LocalFiniteElement, SmallSingletonMatrix> const
frictionalBoundaryMassAssembler(frictionalBoundary);
p1Assembler.assembleOperator(frictionalBoundaryMassAssembler,
frictionalBoundaryMassMatrix);
}
// Q: Does it make sense to weigh them in this manner?
SumNorm<VectorType> const AMNorm(1.0, ANorm, 1.0, MNorm);
SumNorm<Vector> const AMNorm(1.0, ANorm, 1.0, MNorm);
auto const nodalIntegrals =
assembleFrictionWeightsal<GridView, SmallVector>(leafView, p1Assembler,
frictionalNodes);
auto myGlobalNonlinearity = assembleNonlinearity<MatrixType, VectorType>(
auto myGlobalNonlinearity = assembleNonlinearity<Matrix, Vector>(
frictionalNodes, *nodalIntegrals, frictionData);
// Problem formulation: right-hand side
auto const createRHS = [&](double _relativeTime, VectorType &_ell) {
auto const createRHS = [&](double _relativeTime, Vector &_ell) {
assembleNeumann(leafView, p1Assembler, neumannNodes, _ell,
neumannFunction, _relativeTime);
_ell += gravityFunctional;
};
VectorType ell(fineVertexCount);
Vector ell(fineVertexCount);
createRHS(0.0, ell);
// {{{ Initial conditions
SingletonVectorType alpha_initial(fineVertexCount);
SingletonVector alpha_initial(fineVertexCount);
alpha_initial =
std::log(parset.get<double>("boundary.friction.initialState"));
using LinearFactoryType = SolverFactory<
using LinearFactory = SolverFactory<
dims, BlockNonlinearTNNMGProblem<ConvexProblem<
ZeroNonlinearity<SmallVector, SmallMatrix>, MatrixType>>,
GridType>;
ZeroNonlinearity<SmallVector, SmallMatrix>, Matrix>>,
Grid>;
ZeroNonlinearity<SmallVector, SmallMatrix> zeroNonlinearity;
// Solve the stationary problem
VectorType u_initial(fineVertexCount);
Vector u_initial(fineVertexCount);
u_initial = 0.0;
{
LinearFactoryType displacementFactory(parset.sub("solver.tnnmg"), // FIXME
LinearFactory displacementFactory(parset.sub("solver.tnnmg"), // FIXME
refinements, *grid,
displacementDirichletNodes);
auto multigridStep = displacementFactory.getSolver();
typename LinearFactoryType::ConvexProblemType convexProblem(
typename LinearFactory::ConvexProblem convexProblem(
1.0, A, zeroNonlinearity, ell, u_initial);
typename LinearFactoryType::BlockProblemType initialDisplacementProblem(
typename LinearFactory::BlockProblem initialDisplacementProblem(
parset, convexProblem);
auto const &localParset = parset.sub("u0.solver");
multigridStep->setProblem(u_initial, initialDisplacementProblem);
LoopSolver<VectorType> initialDisplacementProblemSolver(
LoopSolver<Vector> initialDisplacementProblemSolver(
multigridStep, localParset.get<size_t>("maximumIterations"),
localParset.get<double>("tolerance"), &ANorm,
localParset.get<Solver::VerbosityMode>("verbosity"),
......@@ -387,7 +378,7 @@ int main(int argc, char *argv[]) {
initialDisplacementProblemSolver.preprocess();
initialDisplacementProblemSolver.solve();
}
VectorType v_initial(fineVertexCount);
Vector v_initial(fineVertexCount);
v_initial = 0.0;
{
// Prescribe a homogeneous velocity field in the x-direction
......@@ -398,13 +389,13 @@ int main(int argc, char *argv[]) {
for (auto &x : v_initial)
x[0] = v_initial_const;
}
VectorType a_initial(fineVertexCount);
Vector a_initial(fineVertexCount);
a_initial = 0.0;
{
/* We solve Au + Cv + Ma + Psi(v) = ell, thus
Ma = - (Au + Cv + Psi(v) - ell)
*/
VectorType accelerationRHS(fineVertexCount);
Vector accelerationRHS(fineVertexCount);
{
accelerationRHS = 0.0;
Arithmetic::addProduct(accelerationRHS, A, u_initial);
......@@ -415,19 +406,19 @@ int main(int argc, char *argv[]) {
accelerationRHS -= ell;
accelerationRHS *= -1.0;
}
LinearFactoryType accelerationFactory(parset.sub("solver.tnnmg"), // FIXME
LinearFactory accelerationFactory(parset.sub("solver.tnnmg"), // FIXME
refinements, *grid,
accelerationDirichletNodes);
auto multigridStep = accelerationFactory.getSolver();
typename LinearFactoryType::ConvexProblemType convexProblem(
typename LinearFactory::ConvexProblem convexProblem(
1.0, M, zeroNonlinearity, accelerationRHS, a_initial);
typename LinearFactoryType::BlockProblemType initialAccelerationProblem(
typename LinearFactory::BlockProblem initialAccelerationProblem(
parset, convexProblem);
auto const &localParset = parset.sub("a0.solver");
multigridStep->setProblem(a_initial, initialAccelerationProblem);
LoopSolver<VectorType> initialAccelerationProblemSolver(
LoopSolver<Vector> initialAccelerationProblemSolver(
multigridStep, localParset.get<size_t>("maximumIterations"),
localParset.get<double>("tolerance"), &MNorm,
localParset.get<Solver::VerbosityMode>("verbosity"),
......@@ -441,12 +432,11 @@ int main(int argc, char *argv[]) {
writer.writeInfo(alpha_initial, u_initial, v_initial);
// Set up TNNMG solver
using NonlinearFactoryType = SolverFactory<
dims,
MyBlockProblem<ConvexProblem<
Dune::GlobalNonlinearity<MatrixType, VectorType>, MatrixType>>,
GridType>;
NonlinearFactoryType factory(parset.sub("solver.tnnmg"), refinements, *grid,
using NonlinearFactory = SolverFactory<
dims, MyBlockProblem<ConvexProblem<
Dune::GlobalNonlinearity<Matrix, Vector>, Matrix>>,
Grid>;
NonlinearFactory factory(parset.sub("solver.tnnmg"), refinements, *grid,
velocityDirichletNodes);
auto multigridStep = factory.getSolver();
......@@ -464,12 +454,12 @@ int main(int argc, char *argv[]) {
initTimeStepper(parset.get<Config::scheme>("timeSteps.scheme"),
velocityDirichletFunction, velocityDirichletNodes, M, A,
C, u_initial, v_initial, a_initial);
auto stateUpdater = initStateUpdater<SingletonVectorType, VectorType>(
auto stateUpdater = initStateUpdater<SingletonVector, Vector>(
parset.get<Config::stateModel>("boundary.friction.stateModel"),
alpha_initial, frictionalNodes, frictionData);
VectorType v = v_initial;
SingletonVectorType alpha(fineVertexCount);
Vector v = v_initial;
SingletonVector alpha(fineVertexCount);
auto const timesteps = parset.get<size_t>("timeSteps.number"),
maximumStateFPI = parset.get<size_t>("v.fpi.maximumIterations"),
......@@ -494,29 +484,29 @@ int main(int argc, char *argv[]) {
auto const relativeTime = double(run) / double(timesteps);
createRHS(relativeTime, ell);
MatrixType velocityMatrix;
VectorType velocityRHS(fineVertexCount);
VectorType velocityIterate(fineVertexCount);
Matrix velocityMatrix;
Vector velocityRHS(fineVertexCount);
Vector velocityIterate(fineVertexCount);
stateUpdater->setup(tau);
timeSteppingScheme->setup(ell, tau, relativeTime, velocityRHS,
velocityIterate, velocityMatrix);
LoopSolver<VectorType> velocityProblemSolver(
multigridStep, maximumIterations, tolerance, &AMNorm, verbosity,
LoopSolver<Vector> velocityProblemSolver(multigridStep, maximumIterations,
tolerance, &AMNorm, verbosity,
false); // absolute error
size_t iterationCounter;
auto solveVelocityProblem = [&](VectorType &_velocityIterate,
SingletonVectorType const &_alpha) {
auto solveVelocityProblem = [&](Vector &_velocityIterate,
SingletonVector const &_alpha) {
myGlobalNonlinearity->updateLogState(_alpha);
// NIT: Do we really need to pass u here?
typename NonlinearFactoryType::ConvexProblemType const convexProblem(
typename NonlinearFactory::ConvexProblem const convexProblem(
1.0, velocityMatrix, *myGlobalNonlinearity, velocityRHS,
_velocityIterate);
typename NonlinearFactoryType::BlockProblemType velocityProblem(
parset, convexProblem);
typename NonlinearFactory::BlockProblem velocityProblem(parset,
convexProblem);
multigridStep->setProblem(_velocityIterate, velocityProblem);
velocityProblemSolver.preprocess();
......@@ -527,9 +517,9 @@ int main(int argc, char *argv[]) {
// Since the velocity explodes in the quasistatic case, use the
// displacement as a convergence criterion
// Q: is this reasonable?
VectorType u;
VectorType u_saved;
SingletonVectorType alpha_saved;
Vector u;
Vector u_saved;
SingletonVector alpha_saved;
double lastStateCorrection;
for (size_t stateFPI = 1; stateFPI <= maximumStateFPI; ++stateFPI) {
stateUpdater->solve(v);
......@@ -580,10 +570,10 @@ int main(int argc, char *argv[]) {
if (parset.get<bool>("io.writeVTK")) {
auto const gridDisplacement =
Dune::make_shared<BasisGridFunction<P1Basis, VectorType> const>(
p1Basis, u);
SingletonVectorType vonMisesStress;
VonMisesStressAssembler<GridType, P0Basis::LocalFiniteElement>
Dune::make_shared<BasisGridFunction<P1Basis, Vector> const>(p1Basis,
u);
SingletonVector vonMisesStress;
VonMisesStressAssembler<Grid, P0Basis::LocalFiniteElement>
localStressAssembler(E, nu, gridDisplacement);
p0Assembler.assembleFunctional(localStressAssembler, vonMisesStress);
......
......@@ -11,9 +11,9 @@
#include "solverfactory.hh"
template <size_t dim, class BlockProblemType, class GridType>
SolverFactory<dim, BlockProblemType, GridType>::SolverFactory(
Dune::ParameterTree const &parset, size_t refinements, GridType const &grid,
template <size_t dim, class BlockProblem, class Grid>
SolverFactory<dim, BlockProblem, Grid>::SolverFactory(
Dune::ParameterTree const &parset, size_t refinements, Grid const &grid,
Dune::BitSetVector<dim> const &ignoreNodes)
: baseEnergyNorm(linearBaseSolverStep),
linearBaseSolver(&linearBaseSolverStep,
......@@ -21,7 +21,7 @@ SolverFactory<dim, BlockProblemType, GridType>::SolverFactory(
parset.get<double>("linear.tolerance"), &baseEnergyNorm,
Solver::QUIET),
transferOperators(refinements),
multigridStep(new SolverType(linearIterationStep, nonlinearSmoother)) {
multigridStep(new Solver(linearIterationStep, nonlinearSmoother)) {
// linear iteration step
linearIterationStep.setMGType(parset.get<int>("linear.cycle"),
parset.get<int>("linear.pre"),
......@@ -31,8 +31,8 @@ SolverFactory<dim, BlockProblemType, GridType>::SolverFactory(
// transfer operators
for (auto &x : transferOperators)
x = new CompressedMultigridTransfer<VectorType>;
TransferOperatorAssembler<GridType>(grid)
x = new CompressedMultigridTransfer<Vector>;
TransferOperatorAssembler<Grid>(grid)
.assembleOperatorPointerHierarchy(transferOperators);
linearIterationStep.setTransferOperators(transferOperators);
......@@ -43,17 +43,16 @@ SolverFactory<dim, BlockProblemType, GridType>::SolverFactory(
multigridStep->ignoreNodes_ = &ignoreNodes;
}
template <size_t dim, class BlockProblemType, class GridType>
SolverFactory<dim, BlockProblemType, GridType>::~SolverFactory() {
template <size_t dim, class BlockProblem, class Grid>
SolverFactory<dim, BlockProblem, Grid>::~SolverFactory() {
for (auto &x : transferOperators)
delete x;
delete multigridStep;
}
template <size_t dim, class BlockProblemType, class GridType>
auto SolverFactory<dim, BlockProblemType, GridType>::getSolver()
-> SolverType *{
template <size_t dim, class BlockProblem, class Grid>
auto SolverFactory<dim, BlockProblem, Grid>::getSolver() -> Solver *{
return multigridStep;
}
......
......@@ -16,38 +16,37 @@
#include <dune/tnnmg/iterationsteps/genericnonlineargs.hh>
#include <dune/tnnmg/iterationsteps/tnnmgstep.hh>
template <size_t dim, class BlockProblemTypeTEMPLATE, class GridType>
template <size_t dim, class BlockProblemTEMPLATE, class Grid>
class SolverFactory {
public:
using BlockProblemType = BlockProblemTypeTEMPLATE;
using ConvexProblemType = typename BlockProblemType::ConvexProblemType;
using VectorType = typename BlockProblemType::VectorType;
using MatrixType = typename BlockProblemType::MatrixType;
using BlockProblem = BlockProblemTEMPLATE;
using ConvexProblem = typename BlockProblem::ConvexProblemType;
using Vector = typename BlockProblem::VectorType;
using Matrix = typename BlockProblem::MatrixType;
private:
using NonlinearSmootherType = GenericNonlinearGS<BlockProblemType>;
using SolverType = TruncatedNonsmoothNewtonMultigrid<BlockProblemType,
NonlinearSmootherType>;
using NonlinearSmoother = GenericNonlinearGS<BlockProblem>;
using Solver =
TruncatedNonsmoothNewtonMultigrid<BlockProblem, NonlinearSmoother>;
public:
SolverFactory(Dune::ParameterTree const &parset, size_t refinements,
GridType const &grid,
Dune::BitSetVector<dim> const &ignoreNodes);
Grid const &grid, Dune::BitSetVector<dim> const &ignoreNodes);
~SolverFactory();
SolverType *getSolver();
Solver *getSolver();
private:
TruncatedBlockGSStep<MatrixType, VectorType> linearBaseSolverStep;
EnergyNorm<MatrixType, VectorType> baseEnergyNorm;
LoopSolver<VectorType> linearBaseSolver;
TruncatedBlockGSStep<MatrixType, VectorType> linearPresmoother;
TruncatedBlockGSStep<MatrixType, VectorType> linearPostsmoother;
MultigridStep<MatrixType, VectorType> linearIterationStep;
std::vector<CompressedMultigridTransfer<VectorType> *> transferOperators;
NonlinearSmootherType nonlinearSmoother;
SolverType *multigridStep;
TruncatedBlockGSStep<Matrix, Vector> linearBaseSolverStep;
EnergyNorm<Matrix, Vector> baseEnergyNorm;
LoopSolver<Vector> linearBaseSolver;
TruncatedBlockGSStep<Matrix, Vector> linearPresmoother;
TruncatedBlockGSStep<Matrix, Vector> linearPostsmoother;
MultigridStep<Matrix, Vector> linearIterationStep;
std::vector<CompressedMultigridTransfer<Vector> *> transferOperators;
NonlinearSmoother nonlinearSmoother;
Solver *multigridStep;
};
#endif
......@@ -22,16 +22,16 @@
using SmallVector = Dune::FieldVector<double, DIM>;
using SmallMatrix = Dune::FieldMatrix<double, DIM, DIM>;
using VectorType = Dune::BlockVector<SmallVector>;
using MatrixType = Dune::BCRSMatrix<SmallMatrix>;
using Vector = Dune::BlockVector<SmallVector>;
using Matrix = Dune::BCRSMatrix<SmallMatrix>;
using GridType = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
using Grid = Dune::ALUGrid<DIM, DIM, Dune::simplex, Dune::nonconforming>;
template class SolverFactory<
DIM, MyBlockProblem<ConvexProblem<
Dune::GlobalNonlinearity<MatrixType, VectorType>, MatrixType>>,
GridType>;
DIM, MyBlockProblem<
ConvexProblem<Dune::GlobalNonlinearity<Matrix, Vector>, Matrix>>,
Grid>;
template class SolverFactory<
DIM, BlockNonlinearTNNMGProblem<ConvexProblem<
ZeroNonlinearity<SmallVector, SmallMatrix>, MatrixType>>,
GridType>;
DIM, BlockNonlinearTNNMGProblem<
ConvexProblem<ZeroNonlinearity<SmallVector, SmallMatrix>, Matrix>>,
Grid>;
......@@ -4,46 +4,44 @@
#include "compute_state_dieterich_euler.hh"
#include "stateupdater.hh"
template <class SingletonVectorType, class VectorType>
class DieterichStateUpdater
: public StateUpdater<SingletonVectorType, VectorType> {
template <class SingletonVector, class Vector>
class DieterichStateUpdater : public StateUpdater<SingletonVector, Vector> {
public:
DieterichStateUpdater(SingletonVectorType _logState_initial,
DieterichStateUpdater(SingletonVector _logState_initial,
Dune::BitSetVector<1> const &_nodes, double _L);
virtual void nextTimeStep();
virtual void setup(double _tau);
virtual void solve(VectorType const &velocity_field);
virtual void extractLogState(SingletonVectorType &);
virtual void solve(Vector const &velocity_field);
virtual void extractLogState(SingletonVector &);
private:
SingletonVectorType logState_o;
SingletonVectorType logState;
SingletonVector logState_o;
SingletonVector logState;
Dune::BitSetVector<1> const &nodes;
double const L;
double tau;
};
template <class SingletonVectorType, class VectorType>
DieterichStateUpdater<SingletonVectorType, VectorType>::DieterichStateUpdater(
SingletonVectorType _logState_initial, Dune::BitSetVector<1> const &_nodes,
template <class SingletonVector, class Vector>
DieterichStateUpdater<SingletonVector, Vector>::DieterichStateUpdater(
SingletonVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
double _L)
: logState(_logState_initial), nodes(_nodes), L(_L) {}
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
template <class SingletonVector, class Vector>
void DieterichStateUpdater<SingletonVector, Vector>::nextTimeStep() {
logState_o = logState;
}
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::setup(
double _tau) {
template <class SingletonVector, class Vector>
void DieterichStateUpdater<SingletonVector, Vector>::setup(double _tau) {
tau = _tau;
}
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::solve(
VectorType const &velocity_field) {
template <class SingletonVector, class Vector>
void DieterichStateUpdater<SingletonVector, Vector>::solve(
Vector const &velocity_field) {
for (size_t i = 0; i < nodes.size(); ++i)
if (nodes[i][0]) {
double const V = velocity_field[i].two_norm();
......@@ -51,9 +49,9 @@ void DieterichStateUpdater<SingletonVectorType, VectorType>::solve(
}
}
template <class SingletonVectorType, class VectorType>
void DieterichStateUpdater<SingletonVectorType, VectorType>::extractLogState(
SingletonVectorType &_logState) {
template <class SingletonVector, class Vector>
void DieterichStateUpdater<SingletonVector, Vector>::extractLogState(
SingletonVector &_logState) {
_logState = logState;
}
......
......@@ -4,44 +4,44 @@
#include "compute_state_ruina.hh"
#include "stateupdater.hh"
template <class SingletonVectorType, class VectorType>
class RuinaStateUpdater : public StateUpdater<SingletonVectorType, VectorType> {
template <class SingletonVector, class Vector>
class RuinaStateUpdater : public StateUpdater<SingletonVector, Vector> {
public:
RuinaStateUpdater(SingletonVectorType _logState_initial,
RuinaStateUpdater(SingletonVector _logState_initial,
Dune::BitSetVector<1> const &_nodes, double _L);
virtual void nextTimeStep();
virtual void setup(double _tau);
virtual void solve(VectorType const &velocity_field);
virtual void extractLogState(SingletonVectorType &);
virtual void solve(Vector const &velocity_field);
virtual void extractLogState(SingletonVector &);
private:
SingletonVectorType logState_o;
SingletonVectorType logState;
SingletonVector logState_o;
SingletonVector logState;
Dune::BitSetVector<1> const &nodes;
double const L;
double tau;
};
template <class SingletonVectorType, class VectorType>
RuinaStateUpdater<SingletonVectorType, VectorType>::RuinaStateUpdater(
SingletonVectorType _logState_initial, Dune::BitSetVector<1> const &_nodes,
template <class SingletonVector, class Vector>
RuinaStateUpdater<SingletonVector, Vector>::RuinaStateUpdater(
SingletonVector _logState_initial, Dune::BitSetVector<1> const &_nodes,
double _L)
: logState(_logState_initial), nodes(_nodes), L(_L) {}
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::nextTimeStep() {
template <class SingletonVector, class Vector>
void RuinaStateUpdater<SingletonVector, Vector>::nextTimeStep() {
logState_o = logState;
}
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::setup(double _tau) {
template <class SingletonVector, class Vector>
void RuinaStateUpdater<SingletonVector, Vector>::setup(double _tau) {
tau = _tau;
}
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::solve(
VectorType const &velocity_field) {
template <class SingletonVector, class Vector>
void RuinaStateUpdater<SingletonVector, Vector>::solve(
Vector const &velocity_field) {
for (size_t i = 0; i < nodes.size(); ++i)
if (nodes[i][0]) {
double const V = velocity_field[i].two_norm();
......@@ -49,9 +49,9 @@ void RuinaStateUpdater<SingletonVectorType, VectorType>::solve(
}
}
template <class SingletonVectorType, class VectorType>
void RuinaStateUpdater<SingletonVectorType, VectorType>::extractLogState(
SingletonVectorType &_logState) {
template <class SingletonVector, class Vector>
void RuinaStateUpdater<SingletonVector, Vector>::extractLogState(
SingletonVector &_logState) {
_logState = logState;
}
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment