Commit 633b910f authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some cleanup

parent 5a874ddf
Pipeline #903 failed with stage
in 2 minutes and 2 seconds
......@@ -8,8 +8,9 @@ DEPENDENCIES=@REQUIRES@
Name: @PACKAGE_NAME@
Version: @VERSION@
Description: amdis module
URL: http://www.github.com/spraetor
Description: AMDiS dune-module
URL: https://gitlab.math.tu-dresden.de/spraetor/dune-amdis
Requires: dune-common dune-geometry dune-localfunctions dune-istl dune-typetree dune-grid dune-functions
Suggests: dune-uggrid dune-alugrid
Libs: -L${libdir}
Cflags: -I${includedir}
......@@ -8,3 +8,4 @@ Version: 0.1
Maintainer: simon.praetorius@tu-dresden.de
#depending on
Depends: dune-common dune-geometry dune-localfunctions dune-istl dune-typetree dune-grid dune-functions
Suggests: dune-uggrid dune-alugrid
......@@ -196,7 +196,7 @@ namespace AMDiS
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
Dune::FieldVector<double,dim> const& quadPos = quad[iq].position();
// The multiplicative factor in the integral transformation formula
const double factor = geometry.integrationElement(quadPos) * quad[iq].weight();
......@@ -243,8 +243,6 @@ namespace AMDiS
for (auto* operatorTerm : zeroOrder)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -286,8 +284,6 @@ namespace AMDiS
for (auto* operatorTerm : firstOrderGrdPhi)
operatorTerm->init(element, quad);
for (size_t iq = 0; iq < quad.size(); ++iq) {
// Position of the current quadrature point in the reference element
const Dune::FieldVector<double,dim>& quadPos = quad[iq].position();
......@@ -517,7 +513,7 @@ namespace AMDiS
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addFOTImpl(Term const& term,
std::size_t i,
size_t i,
FirstOrderType firstOrderType)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, VectorComponent>;
......@@ -542,7 +538,7 @@ namespace AMDiS
template <class MeshView>
template <class Term>
Operator<MeshView>& Operator<MeshView>::addSOTImpl(Term const& term,
std::size_t i, std::size_t j)
size_t i, size_t j)
{
using OpTerm = GenericOperatorTerm<MeshView, Term, MatrixComponent>;
secondOrder.push_back(new OpTerm(term, {i,j}));
......
......@@ -30,8 +30,7 @@ namespace AMDiS
using PointList = std::vector<Dune::QuadraturePoint<double, dim>>;
public:
virtual void init(Element const& element,
PointList const& points) = 0;
virtual void init(Element const& element, PointList const& points) = 0;
virtual double evalZot(size_t iq,
Dune::FieldVector<double,1> const& test,
......@@ -58,7 +57,7 @@ namespace AMDiS
template <class MeshView, class Term, class Traits = tag::none>
class GenericOperatorTerm
: public OperatorTerm<MeshView>
, public OperatorEvaluation
, public OperatorEvaluation // NOTE: better: use static functions
{
using Super = OperatorTerm<MeshView>;
using Element = typename Super::Element;
......@@ -73,8 +72,7 @@ namespace AMDiS
, traits(traits)
{}
virtual void init(Element const& element,
PointList const& points) override
virtual void init(Element const& element, PointList const& points) override
{
term.init(element, points);
......@@ -124,7 +122,7 @@ namespace AMDiS
using value_type = std::decay_t< decltype( std::declval<Term>()[std::declval<size_t>()] ) >;
using _cat = ValueCategory_t<value_type>;
std::vector<value_type> values;
std::vector<value_type> values; // NOTE: maybe caching is not necessary here, since cached already in the term
};
} // end namespace AMDiS
......@@ -26,17 +26,14 @@ namespace AMDiS
public:
/// Constructs a ProblemInstat with prob as its stationary problem.
ProblemInstat(std::string name,
ProblemType& prob)
: ProblemInstatBase(name, NULL),
problemStat(prob)
ProblemInstat(std::string name, ProblemType& prob)
: ProblemInstatBase(name, nullptr)
, problemStat(prob)
{}
ProblemInstat(std::string name,
ProblemType& prob,
ProblemStatBase& initialProb)
: ProblemInstatBase(name, &initialProb),
problemStat(prob)
ProblemInstat(std::string name, ProblemType& prob, ProblemStatBase& initialProb)
: ProblemInstatBase(name, &initialProb)
, problemStat(prob)
{}
/// Initialisation of the problem.
......@@ -63,7 +60,7 @@ namespace AMDiS
/// Returns the I'th component of \ref oldSolution.
template <size_t I = 0>
decltype(auto) getOldSolution(const index_<I> _i = index_<I>())
decltype(auto) getOldSolution(const index_<I> _i = {})
{
return (*oldSolution)[_i];
}
......@@ -76,7 +73,7 @@ namespace AMDiS
ProblemType& problemStat;
/// Solution of the last timestep.
shared_ptr<SystemVectorType> oldSolution;
std::shared_ptr<SystemVectorType> oldSolution;
};
......
......@@ -7,6 +7,8 @@ namespace AMDiS
template <class Traits>
void ProblemInstat<Traits>::transferInitialSolution(AdaptInfo& adaptInfo)
{
AMDIS_FUNCNAME("ProblemInstat::transferInitialSolution()");
test_exit(adaptInfo.getTime() == adaptInfo.getStartTime(),
"after initial solution: time != start time");
problemStat.writeFiles(adaptInfo, true);
......@@ -24,8 +26,6 @@ namespace AMDiS
template <class Traits>
void ProblemInstat<Traits>::initialize(Flag initFlag)
{
AMDIS_FUNCNAME("ProblemInstat::initialize()");
// === create vector for old solution ===
if (initFlag.isSet(INIT_UH_OLD))
createUhOld();
......@@ -35,6 +35,8 @@ namespace AMDiS
template <class Traits>
void ProblemInstat<Traits>::createUhOld()
{
AMDIS_FUNCNAME("ProblemInstat::createUhOld()");
if (oldSolution) {
warning("oldSolution already created\n");
}
......@@ -43,7 +45,7 @@ namespace AMDiS
// create oldSolution
std::vector<std::string> componentNames(size, name + "_uOld");
oldSolution = make_shared<SystemVectorType>(*problemStat.getFeSpaces(), componentNames);
oldSolution = std::make_shared<SystemVectorType>(*problemStat.getFeSpaces(), componentNames);
}
}
......
......@@ -14,8 +14,9 @@ namespace AMDiS
* \brief
* Base class for \ref ProblemInstat.
*/
class ProblemInstatBase : public ProblemTimeInterface,
public ProblemStatBase // NOTE: Why is this derived from ProblemStatBase
class ProblemInstatBase
: public ProblemTimeInterface
, public ProblemStatBase // NOTE: Why is this derived from ProblemStatBase
{
public:
/// Constructor.
......
......@@ -35,7 +35,8 @@ namespace AMDiS
{
template <class Traits>
class ProblemStatSeq : public ProblemStatBase
class ProblemStatSeq
: public ProblemStatBase
{
using Self = ProblemStatSeq;
......@@ -109,15 +110,15 @@ namespace AMDiS
double* factor = NULL,
double* estFactor = NULL);
void addMatrixOperator(shared_ptr<OperatorType> op,
void addMatrixOperator(std::shared_ptr<OperatorType> op,
int i, int j,
double* factor = NULL,
double* estFactor = NULL);
void addMatrixOperator(std::map< std::pair<int,int>,
shared_ptr<OperatorType> > ops);
std::shared_ptr<OperatorType> > ops);
void addMatrixOperator(std::map< std::pair<int,int>,
std::pair<shared_ptr<OperatorType>, bool> > ops);
std::pair<std::shared_ptr<OperatorType>, bool> > ops);
/** @} */
......@@ -128,15 +129,13 @@ namespace AMDiS
double* factor = NULL,
double* estFactor = NULL);
void addVectorOperator(shared_ptr<OperatorType> op,
void addVectorOperator(std::shared_ptr<OperatorType> op,
int i,
double* factor = NULL,
double* estFactor = NULL);
void addVectorOperator(std::map< int,
shared_ptr<OperatorType> > ops);
void addVectorOperator(std::map< int,
std::pair<shared_ptr<OperatorType>, bool> > ops);
void addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops);
void addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops);
/** @} */
/// Adds a Dirichlet boundary condition
......@@ -162,6 +161,7 @@ namespace AMDiS
/// Writes output files.
void writeFiles(AdaptInfo& adaptInfo, bool force = false);
public: // get-methods
/// Returns nr of components \ref nComponents
......@@ -196,9 +196,9 @@ namespace AMDiS
/// Return the \ref linearSolver
auto getSolver() { return linearSolver; }
void setSolver(std::shared_ptr<LinearSolverType> const& solver_)
void setSolver(std::shared_ptr<LinearSolverType> const& solver)
{
linearSolver = solver_;
linearSolver = solver;
}
/// Return the \ref mesh
......@@ -207,7 +207,7 @@ namespace AMDiS
void setMesh(std::shared_ptr<Mesh> const& mesh_)
{
mesh = mesh_;
meshView = make_shared<MeshView>(mesh->leafGridView());
meshView = std::make_shared<MeshView>(mesh->leafGridView());
createFeSpaces();
createMatricesAndVectors();
......@@ -251,7 +251,7 @@ namespace AMDiS
"No mesh name specified for '", name, "->mesh'!");
mesh = MeshCreator<Mesh>::create(meshName);
meshView = make_shared<MeshView>(mesh->leafGridView());
meshView = std::make_shared<MeshView>(mesh->leafGridView());
msg("Create mesh:");
msg("#elements = " , mesh->size(0));
......@@ -262,16 +262,16 @@ namespace AMDiS
void createFeSpaces()
{
feSpaces = make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
feSpaces = std::make_shared<FeSpaces>(constructTuple<FeSpaces>(*meshView));
}
void createMatricesAndVectors()
{
systemMatrix = make_shared<SystemMatrixType>(*feSpaces);
solution = make_shared<SystemVectorType>(*feSpaces, componentNames);
systemMatrix = std::make_shared<SystemMatrixType>(*feSpaces);
solution = std::make_shared<SystemVectorType>(*feSpaces, componentNames);
auto rhsNames = std::vector<std::string>(nComponents, "rhs");
rhs = make_shared<SystemVectorType>(*feSpaces, rhsNames);
rhs = std::make_shared<SystemVectorType>(*feSpaces, rhsNames);
}
void createSolver()
......@@ -287,11 +287,10 @@ namespace AMDiS
void createFileWriter()
{
filewriter = make_shared<FileWriter<Traits>>(name + "->output",
*meshView,
componentNames);
filewriter = std::make_shared<FileWriter<Traits>>(name + "->output", *meshView, componentNames);
}
protected: // sub-methods to assemble DOFMatrix
template <class LhsData, class RhsData, class Elements>
......@@ -301,13 +300,13 @@ namespace AMDiS
bool getElementMatrix(RowView const& rowView,
ColView const& colView,
ElementMatrix& elementMatrix,
std::list<shared_ptr<OperatorType>>& operators,
std::list<std::shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors);
template <class RowView>
bool getElementVector(RowView const& rowView,
ElementVector& elementVector,
std::list<shared_ptr<OperatorType>>& operators,
std::list<std::shared_ptr<OperatorType>>& operators,
std::list<double*> const& factors);
......@@ -322,15 +321,6 @@ namespace AMDiS
IndexSet const& indexSet,
ElementVector const& elementvector);
private: // some internal methods
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.
......@@ -341,35 +331,35 @@ namespace AMDiS
std::vector<std::string> componentNames;
/// Mesh of this problem.
shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
std::shared_ptr<Mesh> mesh; // TODO: generalize to multi-mesh problems
/// Name of the mesh
std::string meshName;
/// A gridView object
shared_ptr<MeshView> meshView;
std::shared_ptr<MeshView> meshView;
/// Pointer to the meshes for the different problem components
std::vector<Mesh*> componentMeshes;
/// FE spaces of this problem.
shared_ptr<FeSpaces> feSpaces; // eventuell const
std::shared_ptr<FeSpaces> feSpaces; // eventuell const
/// A FileWriter object
shared_ptr<FileWriter<Traits>> filewriter;
std::shared_ptr<FileWriter<Traits>> filewriter;
/// An object of the linearSolver Interface
shared_ptr<LinearSolverType> linearSolver;
std::shared_ptr<LinearSolverType> linearSolver;
/// A block-matrix that is filled during assembling
shared_ptr<SystemMatrixType> systemMatrix;
std::shared_ptr<SystemMatrixType> systemMatrix;
/// A block-vector with the solution components
shared_ptr<SystemVectorType> solution;
std::shared_ptr<SystemVectorType> solution;
/// A block-vector (load-vector) corresponding to the right.hand side
/// of the equation, filled during assembling
shared_ptr<SystemVectorType> rhs;
std::shared_ptr<SystemVectorType> rhs;
template <class T>
......@@ -381,48 +371,47 @@ namespace AMDiS
/// A map (i,j) -> list<DirichleBC> string a boundary conditions for
/// each matrix block
MatrixEntries<shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
MatrixEntries<std::shared_ptr<DirichletBC<WorldVector>>> dirichletBc;
/// A map (i,j) -> list<OperatorType> string the differential operators for
/// each matrix block
MatrixEntries<shared_ptr<OperatorType>> matrixOperators;
MatrixEntries<std::shared_ptr<OperatorType>> matrixOperators;
MatrixEntries<double*> matrixFactors;
std::map< std::pair<int,int>, bool > matrixAssembled; // if false, do reassemble
std::map< std::pair<int,int>, bool > matrixChanging; // if true, or vectorAssembled false, do reassemble
/// A map (i) -> list<OperatorType> string the differential operators for
/// each rhs block
VectorEntries<shared_ptr<OperatorType>> vectorOperators;
VectorEntries<std::shared_ptr<OperatorType>> vectorOperators;
VectorEntries<double*> vectorFactors;
std::map< int, bool > vectorAssembled; // if false, do reassemble
std::map< int, bool > vectorChanging; // if true, or vectorAssembled false, do reassemble
template <int r, int c>
template <int R, int C>
struct MatrixData
{
using DOFMatrixType =
std::tuple_element_t<c, std::tuple_element_t<r, typename SystemMatrixType::DOFMatrices>>;
std::tuple_element_t<C, std::tuple_element_t<R, typename SystemMatrixType::DOFMatrices>>;
DOFMatrixType& matrix;
std::list<shared_ptr<OperatorType>>& operators;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble;
std::pair<int,int> row_col = {r, c};
std::pair<int,int> row_col = {R, C};
};
template <int r>
template <int I>
struct VectorData
{
using DOFVectorType =
std::tuple_element_t<r, typename SystemVectorType::DOFVectors>;
using DOFVectorType = std::tuple_element_t<I, typename SystemVectorType::DOFVectors>;
DOFVectorType& vector;
std::list<shared_ptr<OperatorType>>& operators;
std::list<std::shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors;
bool assemble;
int row = r;
int row = I;
};
};
......@@ -437,8 +426,9 @@ namespace AMDiS
namespace Impl
{
template <class ProblemStatType>
struct ProblemStat : public ProblemStatType,
public StandardProblemIteration
struct ProblemStat
: public ProblemStatType
, public StandardProblemIteration
{
using ProblemStatType::getName;
......
......@@ -107,7 +107,7 @@ namespace AMDiS
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(shared_ptr<OperatorType> op,
void ProblemStatSeq<Traits>::addMatrixOperator(std::shared_ptr<OperatorType> op,
int i, int j,
double* factor,
double* estFactor)
......@@ -127,18 +127,18 @@ namespace AMDiS
double* factor,
double* estFactor)
{
addMatrixOperator(make_shared<OperatorType>(op), i, j, factor, estFactor);
addMatrixOperator(std::make_shared<OperatorType>(op), i, j, factor, estFactor);
}
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, shared_ptr<OperatorType> > ops)
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::shared_ptr<OperatorType> > ops)
{
for (auto op : ops)
addMatrixOperator(op.second, op.first.first, op.first.second);
}
template <class Traits>
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<shared_ptr<OperatorType>,bool> > ops)
void ProblemStatSeq<Traits>::addMatrixOperator(std::map< std::pair<int,int>, std::pair<std::shared_ptr<OperatorType>,bool> > ops)
{
for (auto op : ops) {
auto row_col = op.first;
......@@ -149,7 +149,7 @@ namespace AMDiS
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(shared_ptr<OperatorType> op,
void ProblemStatSeq<Traits>::addVectorOperator(std::shared_ptr<OperatorType> op,
int i,
double* factor,
double* estFactor)
......@@ -168,18 +168,18 @@ namespace AMDiS
double* factor,
double* estFactor)
{
addVectorOperator(make_shared<OperatorType>(op), i, factor, estFactor);
addVectorOperator(std::make_shared<OperatorType>(op), i, factor, estFactor);
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, shared_ptr<OperatorType> > ops)
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::shared_ptr<OperatorType> > ops)
{
for (auto op : ops)
addVectorOperator(op.second, op.first);
}
template <class Traits>
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<shared_ptr<OperatorType>, bool> > ops)
void ProblemStatSeq<Traits>::addVectorOperator(std::map< int, std::pair<std::shared_ptr<OperatorType>, bool> > ops)
{
for (auto op : ops) {
vectorOperators[op.first].push_back(op.second.first);
......@@ -259,54 +259,54 @@ namespace AMDiS
Timer t;
For<0, nComponents>::loop([this](auto const _r) {
this->getNodeFactory(_r).initializeIndices();
this->getFeSpace(_r).update(this->getMeshView());
});
size_t nnz = 0;
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const r)
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector_](auto const _r)
{
msg(this->getFeSpace(r).size(), " DOFs for FeSpace[", r, "]");
msg(this->getFeSpace(_r).size(), " DOFs for FeSpace[", _r, "]");
bool asmVector = asmVector_ && (!vectorAssembled[int(r)] || vectorChanging[int(r)]);
bool asmVector = asmVector_ && (!vectorAssembled[int(_r)] || vectorChanging[int(_r)]);
if (asmVector) {
rhs->compress(r);
rhs->getVector(r) = 0.0;
rhs->compress(_r);
rhs->getVector(_r) = 0.0;
}
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector](auto const c)
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector](auto const _c)
{
static const int _c = decltype(c)::value;
static const int _r = decltype(r)::value;
index_<_r> r{};
static const int C = decltype(_c)::value;
static const int R = decltype(_r)::value;
index_<R> _r = {};
using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<_r, _c>;
using VectorData = typename ProblemStatSeq<Traits>::template VectorData<_r>;
using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<R, C>;
using VectorData = typename ProblemStatSeq<Traits>::template VectorData<R>;
// The DOFMatrix which should be assembled
auto& dofmatrix = (*systemMatrix)(r, c);
auto& solution_vec = (*solution)[c];
auto& rhs_vec = (*rhs)[r];
auto& dofmatrix = (*systemMatrix)(_r, _c);
auto& solution_vec = (*solution)[_c];
auto& rhs_vec = (*rhs)[_r];
auto row_col = std::make_pair(_r, _c);
auto row_col = std::make_pair(R, C);
bool assembleMatrix = asmMatrix_ && (!matrixAssembled[row_col] || matrixChanging[row_col]);
bool assembleVector = asmVector && _r == _c;
bool assembleVector = asmVector && R == C;
int r_ = 0, c_ = 0;
if (assembleMatrix) {
// init boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r_, c_) = bc_list.first;
if (r_ == _r) {
if (r_ == R) {
for (auto bc : bc_list.second)
bc->init(c_ == _c, dofmatrix, solution_vec, rhs_vec);
bc->init(c_ == C, dofmatrix, solution_vec, rhs_vec);
}
}
}
auto mat = MatrixData{dofmatrix, matrixOperators[row_col], matrixFactors[row_col], assembleMatrix};
auto vec = VectorData{rhs_vec, vectorOperators[_r], vectorFactors[_r], assembleVector};
auto vec = VectorData{rhs_vec, vectorOperators[R], vectorFactors[R], assembleVector};
// assemble the DOFMatrix block and the corresponding rhs vector, of r==c
......@@ -317,15 +317,15 @@ namespace AMDiS
if (assembleMatrix)
matrixAssembled[row_col] = true;
if (assembleVector)
vectorAssembled[_r] = true;
vectorAssembled[R] = true;
if (assembleMatrix) {
// finish boundary condition
for (auto bc_list : dirichletBc) {
std::tie(r_, c_) = bc_list.first;
if (r_ == _r) {
if (r_ == R) {