Commit 59f0c65c authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

moved to current dune version

parent 9b0f4bf7
Pipeline #901 passed with stage
in 13 minutes and 7 seconds
...@@ -30,8 +30,8 @@ namespace AMDiS ...@@ -30,8 +30,8 @@ namespace AMDiS
{ {
public: public:
template <class Predicate, class Values, template <class Predicate, class Values,
class = std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> && std::enable_if_t< Concepts::Functor<Predicate, bool(WorldVector)> &&
Concepts::Functor<Values, double(WorldVector)> > > Concepts::Functor<Values, double(WorldVector)>, int*> = nullptr>
DirichletBC(Predicate&& predicate, Values&& values) DirichletBC(Predicate&& predicate, Values&& values)
: predicate(std::forward<Predicate>(predicate)) : predicate(std::forward<Predicate>(predicate))
, values(std::forward<Values>(values)) , values(std::forward<Values>(values))
...@@ -40,16 +40,16 @@ namespace AMDiS ...@@ -40,16 +40,16 @@ namespace AMDiS
template <class Matrix, class VectorX, class VectorB> template <class Matrix, class VectorX, class VectorB>
void init(bool apply, void init(bool apply,
Matrix& matrix, Matrix& matrix,
VectorX& rhs, VectorX& rhs,
VectorB& solution); VectorB& solution);
template <class Matrix, class VectorX, class VectorB> template <class Matrix, class VectorX, class VectorB>
void finish(bool apply, void finish(bool apply,
Matrix& matrix, Matrix& matrix,
VectorX& rhs, VectorX& rhs,
VectorB& solution); VectorB& solution);
private: private:
std::function<bool(WorldVector)> predicate; std::function<bool(WorldVector)> predicate;
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
#include <dune/grid/io/file/vtk/vtkwriter.hh> #include <dune/grid/io/file/vtk/vtkwriter.hh>
#include <dune/grid/io/file/vtk/vtksequencewriter.hh> #include <dune/grid/io/file/vtk/vtksequencewriter.hh>
#include <dune/geometry/genericgeometry/referenceelements.hh> #include <dune/geometry/referenceelements.hh>
#include <dune/amdis/Initfile.hpp> #include <dune/amdis/Initfile.hpp>
#include <dune/amdis/common/Loops.hpp> #include <dune/amdis/common/Loops.hpp>
...@@ -19,22 +19,22 @@ namespace AMDiS ...@@ -19,22 +19,22 @@ namespace AMDiS
class FileWriter class FileWriter
{ {
private: // typedefs and static constants private: // typedefs and static constants
using Mesh = typename Traits::Mesh; using Mesh = typename Traits::Mesh;
using MeshView = typename Mesh::LeafGridView; using MeshView = typename Mesh::LeafGridView;
/// Dimension of the mesh /// Dimension of the mesh
static constexpr int dim = Traits::dim; static constexpr int dim = Traits::dim;
/// Dimension of the world /// Dimension of the world
static constexpr int dow = Traits::dimworld; static constexpr int dow = Traits::dimworld;
/// Number of problem components /// Number of problem components
static constexpr int nComponents = Traits::nComponents; static constexpr int nComponents = Traits::nComponents;
public: public:
/// Constructor. /// Constructor.
FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names_) FileWriter(std::string base, MeshView const& meshView, std::vector<std::string> const& names_)
: meshView(meshView) : meshView(meshView)
...@@ -44,26 +44,26 @@ namespace AMDiS ...@@ -44,26 +44,26 @@ namespace AMDiS
Parameters::get(base + "->filename", filename); Parameters::get(base + "->filename", filename);
dir = "output"; dir = "output";
Parameters::get(base + "->output directory", dir); Parameters::get(base + "->output directory", dir);
vtkWriter = make_shared<Dune::VTKWriter<MeshView>>(meshView); vtkWriter = make_shared<Dune::VTKWriter<MeshView>>(meshView);
vtkSeqWriter = make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, ""); vtkSeqWriter = make_shared<Dune::VTKSequenceWriter<MeshView>>(vtkWriter, filename, dir, "");
} }
/// default write method for time-depended data /// default write method for time-depended data
template <class SystemVectorType> template <class SystemVectorType>
void write(double time, SystemVectorType const& solutions) void write(double time, SystemVectorType const& solutions)
{ {
vtkWriter->clear(); vtkWriter->clear();
// copy dofvector to vertex data // copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i) for_<0, nComponents>([this, &solutions](const auto _i)
{ {
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors)); this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]); vtkSeqWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
}); });
vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/); vtkSeqWriter->write(time/*, Dune::VTK::appendedraw*/);
} }
/// default write method for stationary data /// default write method for stationary data
template <class SystemVectorType> template <class SystemVectorType>
...@@ -71,48 +71,48 @@ namespace AMDiS ...@@ -71,48 +71,48 @@ namespace AMDiS
{ {
vtkWriter->clear(); vtkWriter->clear();
// copy dofvector to vertex data // copy dofvector to vertex data
for_<0, nComponents>([this, &solutions](const auto _i) for_<0, nComponents>([this, &solutions](const auto _i)
{ {
this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors)); this->dofVector2vertexVector(solutions[_i], std::get<_i>(data_vectors));
vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]); vtkWriter->addVertexData(std::get<_i>(data_vectors), names[_i]);
}); });
vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/); vtkWriter->pwrite(filename, dir, "" /*, Dune::VTK::appendedraw*/);
} }
protected: protected:
template <class DOFVector, class Vector> template <class DOFVector, class Vector>
void dofVector2vertexVector(DOFVector const& dofvector, Vector& data) void dofVector2vertexVector(DOFVector const& dofvector, Vector& data)
{ {
using Geometry = typename MeshView::template Codim<0>::Geometry; using Geometry = typename MeshView::template Codim<0>::Geometry;
using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>; using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
data.resize(meshView.size(dim)); data.resize(meshView.size(dim));
auto const& indexSet = meshView.indexSet(); auto const& indexSet = meshView.indexSet();
auto const& feSpace = dofvector.getFeSpace(); auto const& feSpace = dofvector.getFeSpace();
auto localView = feSpace.localView(); auto localView = feSpace.localView();
auto localIndexSet = feSpace.localIndexSet(); auto localIndexSet = feSpace.localIndexSet();
// copy data to P1-vector // copy data to P1-vector
for (auto const& element : elements(meshView)) { for (auto const& element : elements(meshView)) {
localView.bind(element); localView.bind(element);
localIndexSet.bind(localView); localIndexSet.bind(localView);
auto const& localBasis = localView.tree().finiteElement().localBasis(); auto const& localBasis = localView.tree().finiteElement().localBasis();
auto const& refElement = RefElements::general(element.type()); auto const& refElement = RefElements::general(element.type());
std::vector<Dune::FieldVector<double,1> > shapeValues; std::vector<Dune::FieldVector<double,1> > shapeValues;
size_t nVertices = element.subEntities(dim); size_t nVertices = element.subEntities(dim);
for (size_t i = 0; i < nVertices; ++i) { for (size_t i = 0; i < nVertices; ++i) {
auto const& v = element.template subEntity<dim>(i); auto const& v = element.template subEntity<dim>(i);
auto pos = refElement.position(i, dim); auto pos = refElement.position(i, dim);
localBasis.evaluateFunction(pos, shapeValues); localBasis.evaluateFunction(pos, shapeValues);
size_t idx = indexSet.index(v); size_t idx = indexSet.index(v);
data[idx] = 0.0; data[idx] = 0.0;
for (size_t j = 0; j < shapeValues.size(); ++j) { for (size_t j = 0; j < shapeValues.size(); ++j) {
...@@ -127,10 +127,10 @@ namespace AMDiS ...@@ -127,10 +127,10 @@ namespace AMDiS
MeshView const& meshView; MeshView const& meshView;
shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter; shared_ptr<Dune::VTKWriter<MeshView>> vtkWriter;
shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter; shared_ptr<Dune::VTKSequenceWriter<MeshView>> vtkSeqWriter;
std::vector<std::string> names; std::vector<std::string> names;
std::array<std::vector<double>, nComponents> data_vectors; std::array<std::vector<double>, nComponents> data_vectors;
std::string filename; std::string filename;
std::string dir; std::string dir;
}; };
......
...@@ -7,7 +7,7 @@ ...@@ -7,7 +7,7 @@
#include <array> #include <array>
#include <memory> #include <memory>
#include <dune/common/array.hh> #include <dune/common/filledarray.hh>
#include <dune/common/fvector.hh> #include <dune/common/fvector.hh>
#include <dune/grid/albertagrid.hh> #include <dune/grid/albertagrid.hh>
...@@ -149,10 +149,10 @@ namespace AMDiS ...@@ -149,10 +149,10 @@ namespace AMDiS
msg("L = ", L); msg("L = ", L);
auto s = Dune::fill_array<int,dim>(2); // number of cells on coarse mesh in each direction auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
Parameters::get(meshName + "->num cells", s); Parameters::get(meshName + "->num cells", s);
// TODO: add more parameters for yasp-grid (see constructor) // TODO: add more parameters for yasp-grid (see constructor)
return make_unique<Grid>(L, s); return make_unique<Grid>(L, s);
} }
...@@ -171,10 +171,10 @@ namespace AMDiS ...@@ -171,10 +171,10 @@ namespace AMDiS
Parameters::get(meshName + "->min corner", lowerleft); Parameters::get(meshName + "->min corner", lowerleft);
Parameters::get(meshName + "->max corner", upperright); Parameters::get(meshName + "->max corner", upperright);
auto s = Dune::fill_array<int,dim>(2); // number of cells on coarse mesh in each direction auto s = Dune::filledArray<std::size_t(dim)>(2); // number of cells on coarse mesh in each direction
Parameters::get(meshName + "->num cells", s); Parameters::get(meshName + "->num cells", s);
// TODO: add more parameters for yasp-grid (see constructor) // TODO: add more parameters for yasp-grid (see constructor)
return make_unique<Grid>(lowerleft, upperright, s); return make_unique<Grid>(lowerleft, upperright, s);
} }
......
...@@ -196,9 +196,25 @@ namespace AMDiS ...@@ -196,9 +196,25 @@ namespace AMDiS
/// Return the \ref linearSolver /// Return the \ref linearSolver
auto getSolver() { return linearSolver; } auto getSolver() { return linearSolver; }
void setSolver(std::shared_ptr<LinearSolverType> const& solver_)
{
linearSolver = solver_;
}
/// Return the \ref mesh /// Return the \ref mesh
auto getMesh() { return mesh; } auto getMesh() { return mesh; }
void setMesh(std::shared_ptr<Mesh> const& mesh_)
{
mesh = mesh_;
meshView = make_shared<MeshView>(mesh->leafGridView());
createFeSpaces();
createMatricesAndVectors();
createFileWriter();
}
/// Return the \ref meshView /// Return the \ref meshView
auto getMeshView() { return meshView; } auto getMeshView() { return meshView; }
...@@ -380,32 +396,32 @@ namespace AMDiS ...@@ -380,32 +396,32 @@ namespace AMDiS
VectorEntries<double*> vectorFactors; VectorEntries<double*> vectorFactors;
std::map< int, bool > vectorAssembled; // if false, do reassemble std::map< int, bool > vectorAssembled; // if false, do reassemble
std::map< int, bool > vectorChanging; // if true, or vectorAssembled 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 struct MatrixData
{ {
using DOFMatrixType = 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; DOFMatrixType& matrix;
std::list<shared_ptr<OperatorType>>& operators; std::list<shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors; std::list<double*> const& factors;
bool assemble; bool assemble;
std::pair<int,int> row_col = {r, c}; std::pair<int,int> row_col = {r, c};
}; };
template <int r> template <int r>
struct VectorData struct VectorData
{ {
using DOFVectorType = using DOFVectorType =
std::tuple_element_t<r, typename SystemVectorType::DOFVectors>; std::tuple_element_t<r, typename SystemVectorType::DOFVectors>;
DOFVectorType& vector; DOFVectorType& vector;
std::list<shared_ptr<OperatorType>>& operators; std::list<shared_ptr<OperatorType>>& operators;
std::list<double*> const& factors; std::list<double*> const& factors;
bool assemble; bool assemble;
int row = r; int row = r;
}; };
}; };
......
...@@ -275,8 +275,9 @@ namespace AMDiS ...@@ -275,8 +275,9 @@ namespace AMDiS
rhs->getVector(_r) = 0.0; rhs->getVector(_r) = 0.0;
} }
For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector, _r](auto const _c) For<0, nComponents>::loop([this, &nnz, asmMatrix_, asmVector, r=_r](auto const _c)
{ {
auto const _r = r;
using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<_r, _c>; using MatrixData = typename ProblemStatSeq<Traits>::template MatrixData<_r, _c>;
using VectorData = typename ProblemStatSeq<Traits>::template VectorData<_r>; using VectorData = typename ProblemStatSeq<Traits>::template VectorData<_r>;
......
...@@ -20,67 +20,67 @@ using StokesProblem = ProblemStat<StokesParam>; ...@@ -20,67 +20,67 @@ using StokesProblem = ProblemStat<StokesParam>;
int main(int argc, char** argv) int main(int argc, char** argv)
{ {
AMDiS::init(argc, argv); AMDiS::init(argc, argv);
StokesProblem prob("stokes"); StokesProblem prob("stokes");
prob.initialize(INIT_ALL); prob.initialize(INIT_ALL);
double viscosity = 1.0; double viscosity = 1.0;
Parameters::get("stokes->viscosity", viscosity); Parameters::get("stokes->viscosity", viscosity);
// define the differential operators // define the differential operators
using Op = StokesProblem::OperatorType; using Op = StokesProblem::OperatorType;
for_<0,DOW>([&](auto _i) for_<0,DOW>([&](auto _i)
{ {
// <viscosity*grad(u_i), grad(v_i)> // <viscosity*grad(u_i), grad(v_i)>
Op* opL = new Op; Op* opL = new Op;
opL->addSOT( constant(viscosity) ); opL->addSOT( constant(viscosity) );
prob.addMatrixOperator(*opL, _i, _i); prob.addMatrixOperator(*opL, _i, _i);
// <p, d_i(v_i)> // <p, d_i(v_i)>
Op* opB = new Op; Op* opB = new Op;
opB->addFOT( constant(1.0), _i, GRD_PSI ); opB->addFOT( constant(1.0), _i, GRD_PSI );
prob.addMatrixOperator(*opB, _i, DOW); prob.addMatrixOperator(*opB, _i, DOW);
// <d_i(u_i), q> // <d_i(u_i), q>
Op* opDiv = new Op; Op* opDiv = new Op;
opDiv->addFOT( constant(1.0), _i, GRD_PHI ); opDiv->addFOT( constant(1.0), _i, GRD_PHI );
prob.addMatrixOperator(*opDiv, DOW, _i); prob.addMatrixOperator(*opDiv, DOW, _i);
}); });
Op* opZero = new Op; Op* opZero = new Op;
opZero->addZOT( constant(0.0) ); opZero->addZOT( constant(0.0) );
prob.addMatrixOperator(*opZero, DOW, DOW); prob.addMatrixOperator(*opZero, DOW, DOW);
// define boundary regions // define boundary regions
auto left = [](auto const& x) { return x[0] < 1.e-8; }; auto left = [](auto const& x) { return x[0] < 1.e-8; };
auto box = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; }; auto box = [](auto const& x) { return x[0] > 1.0 - 1.e-8 || x[1] < 1.e-8 || x[1] > 1.0 - 1.e-8; };
// define boundary values // define boundary values
auto zero = [](auto const& x) { return 0.0; }; auto zero = [](auto const& x) { return 0.0; };
auto parabolic = [](auto const& x) { return x[1]*(1.0 - x[1]); }; auto parabolic = [](auto const& x) { return x[1]*(1.0 - x[1]); };
// set boundary conditions // set boundary conditions
for (size_t i = 0; i < DOW; ++i) for (size_t i = 0; i < DOW; ++i)
prob.addDirichletBC(box, i, i, zero); prob.addDirichletBC(box, i, i, zero);
prob.addDirichletBC(left, 0, 0, zero); prob.addDirichletBC(left, 0, 0, zero);
prob.addDirichletBC(left, 1, 1, parabolic); prob.addDirichletBC(left, 1, 1, parabolic);
prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, 2, 2, zero); prob.addDirichletBC([](auto const& x) { return x[0] < 1.e-8 && x[1] < 1.e-8; }, 2, 2, zero);
// set initial values for the solver (maybe not necessary) // set initial values for the solver (maybe not necessary)
*prob.getSolution() = 0.0; *prob.getSolution() = 0.0;
AdaptInfo adaptInfo("adapt", prob.getNumComponents()); AdaptInfo adaptInfo("adapt", prob.getNumComponents());
// assemble and solve system // assemble and solve system
prob.buildAfterCoarsen(adaptInfo, Flag(0)); prob.buildAfterCoarsen(adaptInfo, Flag(0));
prob.solve(adaptInfo); prob.solve(adaptInfo);
// output solution // output solution
prob.writeFiles(adaptInfo); prob.writeFiles(adaptInfo);
AMDiS::finalize(); AMDiS::finalize();
return 0; return 0;
} }
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment