diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 15799aaa6707c156c45780e88c19aee2cea29826..0ea739be60f460ed62c9c4f8baebdd7938dcde9a 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -31,4 +31,12 @@ add_dependencies(examples if (dune-foamgrid_FOUND) add_amdis_executable(NAME surface.2d SOURCES surface.cc DIM 2 DOW 3 ALBERTA_GRID) add_dependencies(examples surface.2d) -endif () \ No newline at end of file +endif () + +if (NOT HAVE_MTL AND NOT HAVE_EIGEN) + add_amdis_executable(NAME stokes4.2d SOURCES stokes4.cc DIM 2 DOW 2) + add_amdis_executable(NAME pfc1.2d SOURCES pfc1.cc DIM 2 DOW 2) + add_amdis_executable(NAME pfc2.2d SOURCES pfc2.cc DIM 2 DOW 2) + add_amdis_executable(NAME pfc3.2d SOURCES pfc3.cc DIM 2 DOW 2) + add_dependencies(examples stokes4.2d pfc1.2d pfc2.2d pfc3.2d) +endif() diff --git a/examples/PfcPrecon.hpp b/examples/PfcPrecon.hpp new file mode 100644 index 0000000000000000000000000000000000000000..095c830bfedb216099fd20fe16edfab9cb429c05 --- /dev/null +++ b/examples/PfcPrecon.hpp @@ -0,0 +1,188 @@ +#include + +#include +#include + +#include +#include +#include +#include +#include + +namespace AMDiS { + +template +class PfcPrecon + : public Dune::Preconditioner +{ + using Matrix = M; + + // cg solver with preconditioner + template class P> + struct SubSolver + { + using Solver = Dune::CGSolver; + using Operator = Dune::MatrixAdapter; + using Preconditioner = P; + + SubSolver(std::string const& prefix) + { + Parameters::get(prefix + "->max iteration", maxIt_); + Parameters::get(prefix + "->reduction", reduction_); + } + + void update(SubMatrix const& mat) + { + op_ = std::make_unique(mat); + precon_ = std::make_unique(mat, 1, 1.0); + solver_ = std::make_unique(*op_, *precon_, reduction_, maxIt_, 0); + } + + void apply(SubVector& x, SubVector& b) + { + Dune::InverseOperatorResult stat; + solver_->apply(x, b, stat); + } + + private: + int maxIt_ = 10; + double reduction_ = 1.e-3; + std::unique_ptr op_; + std::unique_ptr precon_; + std::unique_ptr solver_; + }; + +public: + PfcPrecon(M const& matrix, double const& tau) + : matrix_(matrix) + , solverM_("pfc->solverM") + , solverMpK1_("pfc->solverMpK1") + , solverMpK2_("pfc->solverMpK2") + , tau_(tau) + {} + + void update() + { + delta_ = std::sqrt(tau_); + + auto const& matrixM = matrix_[0][0]; + auto const& matrixK = matrix_[0][2]; + + solverM_.update(matrixM); + + matrixMpK1_ = matrixM; + matrixMpK1_.axpy(delta_, matrixK); + solverMpK1_.update(matrixMpK1_); + + matrixMpK2_ = matrixM; + matrixMpK2_.axpy(std::sqrt(delta_), matrixK); + solverMpK2_.update(matrixMpK2_); + } + + void pre(X& /*x*/, Y& /*b*/) override { /* do nothing */ } + void post(X& /*x*/) override { /* do nothing */ } + + void apply(X& x, Y const& d) override + { + X y(x); + Y b(d); + + auto const& matrixM = matrix_[0][0]; + auto const& matrixK = matrix_[0][2]; + + // 1. + solverM_.apply(y[0], b[0]); + // 2. + matrixK.usmv(-tau_, y[0], b[1]); + solverMpK1_.apply(y[1], b[1]); + // 3.a + matrixM.mv(y[1], b[0]); + solverMpK2_.apply(x[0], b[0]); + // 3.b + matrixM.mv(x[0], b[1]); + solverMpK2_.apply(x[1], b[1]); + // 4. + x[0] = y[0]; + x[0].axpy(1.0/delta_, y[1]); + x[0].axpy(-1.0/delta_, x[1]); + // 5. + matrixK.usmv(-1.0, x[1], b[2]); + solverM_.apply(x[2], b[2]); + } + + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::sequential; + } + +private: + Matrix const& matrix_; + + using SubMat = TYPEOF(std::declval()[0][0]); + using SubVec = TYPEOF(std::declval()[0]); + + SubMat matrixMpK1_; + SubMat matrixMpK2_; + + SubSolver solverM_; + SubSolver solverMpK1_; + SubSolver solverMpK2_; + + double delta_ = 1.0; + double tau_ = 1.0; +}; + + + +template +struct ISTLPreconCreator, M> + : public ISTLPreconCreatorInterface +{ + using Super = ISTLPreconCreatorInterface; + using Precon = PfcPrecon; + using Self = ISTLPreconCreator; + + struct Creator + : public CreatorInterfaceName> + { + Creator(double const& tau) + : tau_(&tau) + {} + + std::unique_ptr> createWithString(std::string /*prefix*/) override + { + return std::make_unique(*tau_); + } + + double const* tau_; + }; + + ISTLPreconCreator(double const& tau) + : tau_(&tau) + {} + + using PreconBase = Dune::Preconditioner; + std::unique_ptr create(M const& mat) const override + { + auto precon = std::make_unique(mat, *tau_); + precon->update(); + return std::move(precon); + } + + double const* tau_; +}; + + +template +inline void registerPfcPrecon(Basis const& basis, double const& tau) +{ + using B = Blocking_t; + using M = MatrixGenerator_t; + using X = VectorGenerator_t; + + auto pfcPrecon = new typename ISTLPreconCreator, M>::Creator(tau); + CreatorMap>::addCreator("pfc", pfcPrecon); +} + + +} // end namespace AMDiS diff --git a/examples/PfcSolver.hpp b/examples/PfcSolver.hpp new file mode 100644 index 0000000000000000000000000000000000000000..3541786449979b2d8d7414a0d48ff736cb623206 --- /dev/null +++ b/examples/PfcSolver.hpp @@ -0,0 +1,71 @@ +#include + +#include "PfcPrecon.hpp" + +namespace AMDiS +{ + + template + class PfcSolver + : public LinearSolverInterface + { + using Self = PfcSolver; + using Super = LinearSolverInterface; + + using LinOp = Dune::MatrixAdapter; +#if DUNE_VERSION_GT(DUNE_ISTL,2,6) + using Solver = Dune::RestartedFlexibleGMResSolver; +#else + using Solver = Dune::RestartedGMResSolver; +#endif + using Precon = PfcPrecon; + + public: + /// Constructor + explicit PfcSolver(std::string prefix, double const& tau) + : tau_(&tau) + { + Parameters::get(prefix + "->max iteration", maxIter_); + } + + private: + /// Implements \ref LinearSolverInterface::solveSystemImpl() + void solveImpl(Mat const& A, X& x, Y const& b, SolverInfo& solverInfo) override + { + Dune::Timer t; + + if (solverInfo.doCreateMatrixData()) { + linOp_ = std::make_unique(A); + precon_ = std::make_unique(A, *tau_); + solver_ = std::make_unique(*linOp_, *precon_, solverInfo.relTolerance(), 30, maxIter_, solverInfo.info()); + } + + precon_->update(); + if (solverInfo.info() > 0) + msg("creation of istl solver needed {} seconds", t.elapsed()); + + Dune::InverseOperatorResult statistics; + + Y b_(b); + solver_->apply(x, b_, statistics); + solverInfo.setError(statistics.converged ? 0 : 1); + solverInfo.setRelResidual(statistics.reduction); + + if (!solverInfo.doStoreMatrixData()) { + linOp_.reset(); + precon_.reset(); + solver_.reset(); + } + + } + + private: + double const* tau_; + int maxIter_ = 1000; + + std::unique_ptr linOp_; + std::unique_ptr precon_; + std::unique_ptr solver_; + }; + +} // end namespace AMDiS \ No newline at end of file diff --git a/examples/init/pfc.dat.2d b/examples/init/pfc.dat.2d new file mode 100644 index 0000000000000000000000000000000000000000..18571bf534326fe7a7ed2719c3c95851755a9e3b --- /dev/null +++ b/examples/init/pfc.dat.2d @@ -0,0 +1,27 @@ +pfcMesh->global refinements: 4 +pfcMesh->num cells: 4 4 +pfcMesh->dimension: 200 200 + +pfc->mesh: pfcMesh + +pfc->solver->name: fgmres +pfc->solver->max iteration: 10000 +pfc->solver->relative tolerance: 1.e-10 +pfc->solver->info: 2 +pfc->solver->precon->name: pfc + +pfc->solver1->max iteration: 10 +pfc->solver1->reduction: 1.e-5 + +pfc->solver2->max iteration: 10 +pfc->solver2->reduction: 1.e-5 + +pfc->output[1]->filename: pfc.2d +pfc->output[1]->output directory: output +pfc->output[1]->ParaView format: 1 +pfc->output[1]->ParaView animation: 1 +pfc->output[1]->ParaView mode: 1 + +adapt->timestep: 0.001 +adapt->start time: 0.0 +adapt->end time: 10.0 diff --git a/examples/init/stokes.dat.2d b/examples/init/stokes.dat.2d index 54542825dba336eba2af96e3b3371f548ff8bad6..53d92c96df3704464ea1fd53c7929461c6d2b137 100644 --- a/examples/init/stokes.dat.2d +++ b/examples/init/stokes.dat.2d @@ -1,10 +1,10 @@ -stokesMesh->global refinements: 0 +stokesMesh->global refinements: 4 stokesMesh->dimension: 1.0 1.0 stokesMesh->num cells: 4 4 stokes->mesh: stokesMesh -stokes->solver->name: bicgstag +stokes->solver->name: bicgstab stokes->solver->ell: 3 stokes->solver->max iteration: 1000 stokes->solver->relative tolerance: 1e-8 @@ -12,6 +12,8 @@ stokes->solver->restart: 50 stokes->solver->precon->name: ilu stokes->solver->info: 2 +stokes->innersolver->max iteration: 2 + stokes->output[0]->filename: stokes_u.2d stokes->output[0]->name: u stokes->output[0]->subsampling: 2 diff --git a/examples/pfc1.cc b/examples/pfc1.cc new file mode 100644 index 0000000000000000000000000000000000000000..0d3b7ff2244acead4071e5b9b01c8a756889ed74 --- /dev/null +++ b/examples/pfc1.cc @@ -0,0 +1,110 @@ +#include +#include +#include +#include + +#include "PfcPrecon.hpp" + +using namespace AMDiS; + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + auto grid = MeshCreator>::create("pfcMesh"); + + using namespace Dune::Functions::BasisFactory; + auto basis = makeBasis(grid->leafGridView(), power<3>(lagrange<1>(), blockedLexicographic()) ); + + auto prob = makeProblemStat("pfc", *grid, basis); + prob.initialize(INIT_ALL); + + AdaptInfo adaptInfo("adapt"); + + double r = -0.4; + Parameters::get("pfc->r", r); + + double psi_mean = -0.3; + Parameters::get("pfc->psi_mean", psi_mean); + + auto psi = prob.solution(1); + auto M = 1.0; //max(1.e-5, psi + 1.5); + auto tau = std::ref(adaptInfo.timestep()); + + // mu + prob.addMatrixOperator(zot(1.0), 0, 0); + // -(1+r)*psi + prob.addMatrixOperator(zot(-(1.0+r)), 0, 1); + // -3*psi'^2*psi + prob.addMatrixOperator(zot(-3.0*pow<2>(psi)), 0, 1); + // -2*laplace(psi) + prob.addMatrixOperator(sot(2.0), 0, 1); + // -laplace(nu) + prob.addMatrixOperator(sot(1.0), 0, 2); + // -2*psi'^3 + prob.addVectorOperator(zot(-2.0*pow<3>(psi)), 0); + + // psi + prob.addMatrixOperator(zot(1.0), 1, 1); + // -tau*div(M(psi')grad(mu)) + prob.addMatrixOperator(sot(tau*M,1), 1, 0); + // psi' + prob.addVectorOperator(zot(psi), 1); + + // -laplace(psi) + prob.addMatrixOperator(sot(1.0), 2, 1); + // nu + prob.addMatrixOperator(zot(1.0), 2, 2); + + // create the linear solver + // @{ + + using MatrixType = TYPEOF(prob.systemMatrix()->matrix()); + using VectorType = TYPEOF(prob.solutionVector()->vector()); + + // Turn the matrix into a linear operator + Dune::MatrixAdapter linOp(prob.systemMatrix()->matrix()); + + // Define a pfc block preconditioner + PfcPrecon precon(prob.systemMatrix()->matrix(), adaptInfo.timestep()); + + int solverInfo = 1; + Parameters::get("pfc->solver->info", solverInfo); + + // Construct the actual iterative solver +#if DUNE_VERSION_GT(DUNE_ISTL,2,6) + Dune::RestartedFlexibleGMResSolver solver(linOp, precon, 1e-7, 30, 10000, solverInfo); +#else + Dune::RestartedGMResSolver solver(linOp, precon, 1e-7, 30, 10000, solverInfo); +#endif + + // @} + + + std::srand(12345); + psi << [psi_mean](auto const& x) { + return psi_mean + 0.2*(double(std::rand())/RAND_MAX - 0.5); + }; + + prob.writeFiles(adaptInfo); + while (!adaptInfo.reachedEndTime()) + { + // update time + adaptInfo.setTime(adaptInfo.time() + adaptInfo.timestep()); + + // assemble the linear system + prob.assemble(adaptInfo); + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult statistics; + + // Solve the linear system + precon.update(); + solver.apply(prob.solutionVector()->vector(), prob.rhsVector()->vector(), statistics); + + // output solution + prob.writeFiles(adaptInfo); + } + + return 0; +} diff --git a/examples/pfc2.cc b/examples/pfc2.cc new file mode 100644 index 0000000000000000000000000000000000000000..9ac731794b44d10f1fccaa19fa00a7d82bf26dd7 --- /dev/null +++ b/examples/pfc2.cc @@ -0,0 +1,80 @@ +#include +#include +#include +#include +#include +#include + +#include "PfcPrecon.hpp" + +using namespace AMDiS; + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + auto grid = MeshCreator>::create("pfcMesh"); + + using namespace Dune::Functions::BasisFactory; + auto basis = makeBasis(grid->leafGridView(), power<3>(lagrange<1>(), blockedLexicographic()) ); + + AdaptInfo adaptInfo("adapt"); + + // Register the linear solver so that it can be selected in the initfile. + // This must be done before prob.initialize() is called. + registerPfcPrecon(basis,adaptInfo.timestep()); + + auto prob = makeProblemStat("pfc", *grid, basis); + prob.initialize(INIT_ALL); + + auto probInstat = makeProblemInstat("pfc", prob); + probInstat.initialize(INIT_UH_OLD); + + double r = -0.4; + Parameters::get("pfc->r", r); + + double psi_mean = -0.3; + Parameters::get("pfc->psi_mean", psi_mean); + + auto psi = prob.solution(1); + auto psiOld = probInstat.oldSolution(1); + + auto M = 1.0; //max(1.e-5, psi + 1.5); + auto tau = std::ref(adaptInfo.timestep()); + + // mu + prob.addMatrixOperator(zot(1.0), 0, 0); + // -(1+r)*psi + prob.addMatrixOperator(zot(-(1.0+r)), 0, 1); + // -3*psi'^2*psi + prob.addMatrixOperator(zot(-3.0*pow<2>(psi)), 0, 1); + // -2*laplace(psi) + prob.addMatrixOperator(sot(2.0), 0, 1); + // -laplace(nu) + prob.addMatrixOperator(sot(1.0), 0, 2); + // -2*psi'^3 + prob.addVectorOperator(zot(-2.0*pow<3>(psi)), 0); + + // psi + prob.addMatrixOperator(zot(1.0), 1, 1); + // -tau*div(M(psi')grad(mu)) + prob.addMatrixOperator(sot(tau*M,1), 1, 0); + // psi' + prob.addVectorOperator(zot(psi), 1); + + // -laplace(psi) + prob.addMatrixOperator(sot(1.0), 2, 1); + // nu + prob.addMatrixOperator(zot(1.0), 2, 2); + + // set initial condition + std::srand(12345); + psi << [psi_mean](auto const& x) { + return psi_mean + 0.2*(double(std::rand())/RAND_MAX - 0.5); + }; + + AdaptInstationary adaptInstat("adapt", prob, adaptInfo, probInstat, adaptInfo); + adaptInstat.adapt(); + + return 0; +} diff --git a/examples/pfc3.cc b/examples/pfc3.cc new file mode 100644 index 0000000000000000000000000000000000000000..3cfb4399418bd199bcf0c77e775034459a155573 --- /dev/null +++ b/examples/pfc3.cc @@ -0,0 +1,82 @@ +#include +#include +#include +#include +#include +#include + +#include "PfcSolver.hpp" + +using namespace AMDiS; + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + auto grid = MeshCreator>::create("pfcMesh"); + + using namespace Dune::Functions::BasisFactory; + auto basis = makeBasis(grid->leafGridView(), power<3>(lagrange<1>(), blockedLexicographic()) ); + + AdaptInfo adaptInfo("adapt"); + + auto prob = makeProblemStat("pfc", *grid, basis); + prob.initialize(INIT_ALL); + + // register a solver in the problem after initialization + using MatrixType = TYPEOF(prob.systemMatrix()->matrix()); + using VectorType = TYPEOF(prob.solutionVector()->vector()); + PfcSolver solver("pfc->solver", adaptInfo.timestep()); + prob.setSolver(solver); + + auto probInstat = makeProblemInstat("pfc", prob); + probInstat.initialize(INIT_UH_OLD); + + double r = -0.4; + Parameters::get("pfc->r", r); + + double psi_mean = -0.3; + Parameters::get("pfc->psi_mean", psi_mean); + + auto psi = prob.solution(1); + auto psiOld = probInstat.oldSolution(1); + + auto M = 1.0; //max(1.e-5, psi + 1.5); + auto tau = std::ref(adaptInfo.timestep()); + + // mu + prob.addMatrixOperator(zot(1.0), 0, 0); + // -(1+r)*psi + prob.addMatrixOperator(zot(-(1.0+r)), 0, 1); + // -3*psi'^2*psi + prob.addMatrixOperator(zot(-3.0*pow<2>(psi)), 0, 1); + // -2*laplace(psi) + prob.addMatrixOperator(sot(2.0), 0, 1); + // -laplace(nu) + prob.addMatrixOperator(sot(1.0), 0, 2); + // -2*psi'^3 + prob.addVectorOperator(zot(-2.0*pow<3>(psi)), 0); + + // psi + prob.addMatrixOperator(zot(1.0), 1, 1); + // -tau*div(M(psi')grad(mu)) + prob.addMatrixOperator(sot(tau*M,1), 1, 0); + // psi' + prob.addVectorOperator(zot(psi), 1); + + // -laplace(psi) + prob.addMatrixOperator(sot(1.0), 2, 1); + // nu + prob.addMatrixOperator(zot(1.0), 2, 2); + + // set initial condition + std::srand(12345); + psi << [psi_mean](auto const& x) { + return psi_mean + 0.2*(double(std::rand())/RAND_MAX - 0.5); + }; + + AdaptInstationary adaptInstat("adapt", prob, adaptInfo, probInstat, adaptInfo); + adaptInstat.adapt(); + + return 0; +} diff --git a/examples/stokes4.cc b/examples/stokes4.cc new file mode 100644 index 0000000000000000000000000000000000000000..f5862ee39fa4eb6680693401393464bd2a42d8b3 --- /dev/null +++ b/examples/stokes4.cc @@ -0,0 +1,167 @@ +#include +#include +#include +#include + +#include +#include + +using namespace AMDiS; + +template +class StokesPrecon + : public Dune::Preconditioner +{ + using Matrix = M; + + template + struct SubSolver + { + using Solver = Dune::GeneralizedPCGSolver; + + using Operator = Dune::MatrixAdapter; + using CriterionBase = Dune::Amg::AggregationCriterion >; + using Criterion = Dune::Amg::CoarsenCriterion; + using Smoother = Dune::SeqSSOR; + using SmootherArgs = typename Dune::Amg::SmootherTraits::Arguments; + using Preconditioner = Dune::Amg::FastAMG; + + SubSolver(SubMatrix const& mat) + : op_(mat) + { + int maxIt = 1; + Parameters::get("stokes->innersolver->max iteration", maxIt); + + criterion_.setDebugLevel(0); + parms_.setNoPreSmoothSteps(1); + parms_.setNoPostSmoothSteps(1); + parms_.setDebugLevel(0); + precon_ = std::make_unique(op_, criterion_, parms_); + solver_ = std::make_unique(op_, *precon_, 1.e-3, maxIt, 0); + } + + void apply(SubVector& x, SubVector const& b) + { + Dune::InverseOperatorResult stat; + SubVector b_(b); + solver_->apply(x, b_, stat); + } + + Operator op_; + Criterion criterion_; + SmootherArgs smootherArgs_; + Dune::Amg::Parameters parms_; + std::unique_ptr precon_; + std::unique_ptr solver_; + }; + + +public: + template + StokesPrecon(M const& matrix, Q const& pressureMassMatrix) + : velSolver_(matrix[zero_][zero_][0][0]) + , preSolver_(pressureMassMatrix) + {} + + void pre(X& /*x*/, Y& /*b*/) override { /* do nothing */ } + void post(X& /*x*/) override { /* do nothing */ } + + void apply(X& v, Y const& d) override + { + for (int i = 0; i < 2; ++i) + velSolver_.apply(v[zero_][i], d[zero_][i]); + preSolver_.apply(v[one_], d[one_]); + } + + Dune::SolverCategory::Category category() const override + { + return Dune::SolverCategory::sequential; + } + +private: + using VelMatrix = TYPEOF(std::declval()[zero_][zero_][0][0]); + using VelVector = TYPEOF(std::declval()[zero_][0]); + SubSolver velSolver_; + + using PreMatrix = TYPEOF(std::declval()[one_][one_]); + using PreVector = TYPEOF(std::declval()[one_]); + SubSolver preSolver_; +}; + + + +int main(int argc, char** argv) +{ + Environment env(argc, argv); + + using namespace Dune::Functions::BasisFactory; + + Dune::YaspGrid<2> grid({1.0, 1.0}, {4,4}); + auto basis = makeBasis(grid.leafGridView(), + composite(power<2>(lagrange<2>(), blockedLexicographic()), lagrange<1>(), blockedLexicographic()) ); + + auto prob = makeProblemStat("stokes", grid, basis); + prob.initialize(INIT_ALL); + + double viscosity = 1.0; + Parameters::get("stokes->viscosity", viscosity); + + // tree-paths for components + auto _v = Dune::Indices::_0; + + auto opStokes = makeOperator(tag::stokes{}, viscosity); + prob.addMatrixOperator(opStokes); + + // define boundary values + auto parabolic_y = [](auto const& x) + { + return FieldVector{0.0, x[1]*(1.0 - x[1])}; + }; + + FieldVector zero(0); + + // set boundary conditions for velocity + prob.boundaryManager()->setBoxBoundary({1,2,2,2}); + prob.addDirichletBC(1, _v, _v, parabolic_y); + prob.addDirichletBC(2, _v, _v, zero); + + AdaptInfo adaptInfo("adapt"); + + // assemble and solve system + prob.assemble(adaptInfo); + + // assemble a mass matrix on the pressure space + auto pBasis = makeBasis(grid.leafGridView(), lagrange<1>()); + DOFMatrix pMassMatrix(pBasis, pBasis); + using Element = TYPEOF(grid)::Codim<0>::Entity; + pMassMatrix.addOperator(tag::element_operator{}, zot(1.0)); + pMassMatrix.assemble(); + + + // create and apply an ISTL solver + using MatrixType = TYPEOF(prob.systemMatrix()->matrix()); + using VectorType = TYPEOF(prob.solutionVector()->vector()); + + // Technicality: turn the matrix into a linear operator + Dune::MatrixAdapter linOp(prob.systemMatrix()->matrix()); + StokesPrecon precon(prob.systemMatrix()->matrix(), pMassMatrix.matrix()); + + // Construct the actual iterative solver + Dune::MINRESSolver solver( + linOp, // operator to invert + precon, // preconditioner for interation + 1e-10, // desired residual reduction factor + 500, // maximum number of iterations + 2); // verbosity of the solver + + // Object storing some statistics about the solving process + Dune::InverseOperatorResult statistics; + + // Solve the linear system + solver.apply(prob.solutionVector()->vector(), prob.rhsVector()->vector(), statistics); + + // output solution + prob.writeFiles(adaptInfo); + + return 0; +} diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp index c097c9bd4fa21d6269e0d7c2b1800119c0445739..1a7a1efc5e9e4affd520a05e5a41bdbac411566e 100644 --- a/src/amdis/DirichletBC.hpp +++ b/src/amdis/DirichletBC.hpp @@ -3,6 +3,7 @@ #include #include #include +#include #include #include @@ -10,6 +11,8 @@ #include #include #include +#include +#include #include #include #include @@ -31,26 +34,32 @@ namespace AMDiS * support this symmetric modification. As a result, this method returns a list * of columns values, that should be subtracted from the rhs. **/ - template + template class DirichletBC : public BoundaryCondition { using Super = BoundaryCondition; + using MultiIndex = MI; + + template + using BitVector = VectorGenerator_t; public: template ) > - DirichletBC(BM&& boundaryManager, BoundaryType id, Values&& values) + DirichletBC(BM&& boundaryManager, BoundaryType id, Values&& values, bool symmetric = false) : Super(FWD(boundaryManager), id) , values_(FWD(values)) + , symmetric_(symmetric) {} template ), REQUIRES(Concepts::Functor)> - DirichletBC(Predicate&& predicate, Values&& values) + DirichletBC(Predicate&& predicate, Values&& values, bool symmetric = false) : predicate_(FWD(predicate)) , values_(FWD(values)) + , symmetric_(symmetric) {} template @@ -61,8 +70,8 @@ namespace AMDiS } // fill \ref dirichletNodes_ with 1 or 0 whether DOF corresponds to the boundary or not. - template - void init(RB const& rowBasis, CB const& colBasis); + template + void init(RB const& rowBasis, CB const& colBasis, MatIndexSet& matIndexSet); /// \brief Apply dirichlet BC to matrix and vector, i.e., add a unit-row /// to the matrix and optionally delete the corresponding matrix-column. @@ -79,7 +88,9 @@ namespace AMDiS private: Dune::Std::optional> predicate_; std::function values_; - std::vector dirichletNodes_; + std::vector dirichletNodes_; + + bool symmetric_ = false; }; @@ -89,7 +100,7 @@ namespace AMDiS using WorldVector = typename Basis::GridView::template Codim<0>::Geometry::GlobalCoordinate; template - using type = std::list< std::shared_ptr>> >; + using type = std::list< std::shared_ptr, typename Basis::MultiIndex>> >; }; template diff --git a/src/amdis/DirichletBC.inc.hpp b/src/amdis/DirichletBC.inc.hpp index 6f285d24e0cb7c5293d335b1330830bca886b86f..666cdceb28a88f113c37bb4d636ef76fb21a7f36 100644 --- a/src/amdis/DirichletBC.inc.hpp +++ b/src/amdis/DirichletBC.inc.hpp @@ -8,45 +8,87 @@ #include #include +#include namespace AMDiS { -template - template -void DirichletBC:: -init(RB const& rowBasis, CB const& colBasis) +template + template +void DirichletBC:: +init(RB const& rowBasis, CB const& colBasis, MatIndexSet& matIndexSet) { using Dune::Functions::forEachBoundaryDOF; using LV = typename CB::LocalView; using IS = typename CB::GridView::Intersection; - dirichletNodes_.resize(colBasis.dimension()); + + BitVector> bitVector; + vectorWrapper(bitVector).resize(rowBasis); + + dirichletNodes_.clear(); forEachBoundaryDOF(colBasis, [&](int localIndex, LV const& localView, IS const& intersection) { - dirichletNodes_[localView.index(localIndex)] = onBoundary(intersection); + auto idx = localView.index(localIndex); + auto& visited = Dune::Functions::resolveDynamicMultiIndex(bitVector,idx); + if (onBoundary(intersection) && !visited) { + dirichletNodes_.push_back(idx); + visited = 1; + } }); + + // insert diagonal entry into non-zero pattern + for (auto const& index : dirichletNodes_) + matIndexSet.insert(index,index,1); } -template +template template -void DirichletBC:: +void DirichletBC:: fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath) { - auto columns = Constraints::dirichletBC(matrix, dirichletNodes_); + std::vector> columns; - Dune::Hybrid::ifElse(std::is_same, R>{}, - [&](auto id) { - auto subBasis = Dune::Functions::subspaceBasis(*rhs.basis(), rowTreePath); - Dune::Functions::interpolate(subBasis, rhs, values_, dirichletNodes_); + // clear row (and column) of dirichlet nodes + if (symmetric_) { + FakeContainer col; + for (auto const& index : dirichletNodes_) { + MultiIndex row{}; + matrixWrapper(matrix.matrix()).clear_col(row, index, columns); + matrixWrapper(matrix.matrix()).clear_row(index, col); + } + } else { + FakeContainer col; + for (auto const& index : dirichletNodes_) + matrixWrapper(matrix.matrix()).clear_row(index, col); + } + + // insert diagonal 1 value + { auto inserter = matrixInserter(matrix.matrix()); + for (auto const& index : dirichletNodes_) + inserter.insert(index, index, 1); + } + + // create indicator for interpolation + BitVector bitVector; + vectorWrapper(bitVector).resize(*matrix.rowBasis()); + for (auto const& node : dirichletNodes_) + Dune::Functions::resolveDynamicMultiIndex(bitVector,node) = 1; + + // set dirichlet values in solution vector + Dune::Hybrid::ifElse(std::is_same, R>{}, + [&](auto id) -> void { + auto subBasis = Dune::Functions::subspaceBasis(*solution.basis(), colTreePath); + Dune::Functions::interpolate(subBasis, solution, values_, bitVector); }); // subtract columns of dirichlet nodes from rhs for (auto const& triplet : columns) rhs[triplet.row] -= triplet.value * solution[triplet.col]; - Dune::Hybrid::ifElse(std::is_same, R>{}, - [&](auto id) { - auto subBasis = Dune::Functions::subspaceBasis(*solution.basis(), colTreePath); - Dune::Functions::interpolate(subBasis, solution, values_, dirichletNodes_); + // set dirichlet values in rhs vector + Dune::Hybrid::ifElse(std::is_same, R>{}, + [&](auto id) -> void { + auto subBasis = Dune::Functions::subspaceBasis(*rhs.basis(), rowTreePath); + Dune::Functions::interpolate(subBasis, rhs, values_, bitVector); }); } diff --git a/src/amdis/LinearAlgebra.hpp b/src/amdis/LinearAlgebra.hpp index 9a32bd9a09865918d3e0f4b569ea598500eb183c..b7e725a3656a143046a34c6a827fd0c14a113e06 100644 --- a/src/amdis/LinearAlgebra.hpp +++ b/src/amdis/LinearAlgebra.hpp @@ -5,7 +5,6 @@ #if HAVE_MTL -#include #include #include #include @@ -13,14 +12,12 @@ #elif HAVE_EIGEN -#include #include #include #include #else // ISTL -#include #include #include #include diff --git a/src/amdis/Output.hpp b/src/amdis/Output.hpp index 91264e57eb44d532395bb0a7cfb227f27d408c07..b65dd123754cc594e63b0325770d9d4d6211f026 100644 --- a/src/amdis/Output.hpp +++ b/src/amdis/Output.hpp @@ -195,4 +195,10 @@ namespace AMDiS void test_exit_dbg(bool, Args&&...) {} #endif + template + void print_type(std::string prefix = "") + { + std::cout << prefix << __PRETTY_FUNCTION__ << "\n"; + } + } // end namespace AMDiS diff --git a/src/amdis/PeriodicBC.hpp b/src/amdis/PeriodicBC.hpp index c4253a2f9b34b861a7c35546f5f22bbc590d0871..62e969216ddc3ddd564d8e347db06a4f7c0df9bd 100644 --- a/src/amdis/PeriodicBC.hpp +++ b/src/amdis/PeriodicBC.hpp @@ -9,7 +9,8 @@ #include #include #include - +#include +#include namespace AMDiS { @@ -50,6 +51,9 @@ namespace AMDiS { using Super = BoundaryCondition; + template + using BitVector = VectorGenerator_t; + public: using Domain = D; using MultiIndex = MI; @@ -62,8 +66,8 @@ namespace AMDiS , faceTrafo_(std::move(faceTrafo)) {} - template - void init(RB const& rowBasis, CB const& colBasis); + template + void init(RB const& rowBasis, CB const& colBasis, MatIndexSet& matIndexSet); template void fillBoundaryCondition(Mat& A, Sol& x, Rhs& b, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath); @@ -92,9 +96,11 @@ namespace AMDiS private: FaceTrafo faceTrafo_; - std::vector periodicNodes_; - std::map associations_; + std::vector periodicNodes_; + std::map associations_; Dune::Std::optional hasConnectedIntersections_; + + bool symmetric_ = true; }; diff --git a/src/amdis/PeriodicBC.inc.hpp b/src/amdis/PeriodicBC.inc.hpp index b8901f9f1bfdb5f65fa00b70907bbe6af80c0b09..1fe0d4d34bb8d40ab71cf00371c435101c222326 100644 --- a/src/amdis/PeriodicBC.inc.hpp +++ b/src/amdis/PeriodicBC.inc.hpp @@ -6,14 +6,17 @@ #include #include #include + #include +#include +#include namespace AMDiS { template - template + template void PeriodicBC:: -init(RB const& basis, CB const& /*colBasis*/) +init(RB const& basis, CB const& /*colBasis*/, MatIndexSet& /*matIndexSet*/) { if (!bool(hasConnectedIntersections_)) { hasConnectedIntersections_ = true; @@ -35,6 +38,8 @@ init(RB const& basis, CB const& /*colBasis*/) initAssociations(basis); else initAssociations2(basis); + + // TODO: fill matrixIndexSet with associated indices } @@ -46,7 +51,7 @@ initAssociations(Basis const& basis) using std::sqrt; static const auto tol = sqrt(std::numeric_limits::epsilon()); periodicNodes_.clear(); - periodicNodes_.resize(basis.dimension(), false); + //periodicNodes_.resize(basis.dimension(), false); auto localView = basis.localView(); auto seDOFs = subEntityDOFs(basis); @@ -85,7 +90,7 @@ initAssociations(Basis const& basis) for (std::size_t j = 0; j < outsideCoords.size(); ++j) { auto const& y = outsideCoords[j]; if (distance(x,y) < tol) { - periodicNodes_[insideGlobalDOFs[i]] = true; + periodicNodes_.push_back(insideGlobalDOFs[i]); associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]); } } @@ -132,7 +137,7 @@ initAssociations2(Basis const& basis) using std::sqrt; static const auto tol = sqrt(std::numeric_limits::epsilon()); periodicNodes_.clear(); - periodicNodes_.resize(basis.dimension(), false); + //periodicNodes_.resize(basis.dimension(), false); // Dune::ReservedVector,2> using EntitySeed = typename Basis::GridView::Grid::template Codim<0>::EntitySeed; @@ -184,7 +189,7 @@ initAssociations2(Basis const& basis) for (std::size_t j = 0; j < outsideCoords.size(); ++j) { auto const& y = outsideCoords[j]; if (distance(x,y) < tol) { - periodicNodes_[insideGlobalDOFs[i]] = true; + periodicNodes_.push_back(insideGlobalDOFs[i]); associations_[insideGlobalDOFs[i]] = localView.index(outsideDOFs[j]); } } @@ -238,13 +243,69 @@ template void PeriodicBC:: fillBoundaryCondition(Mat& matrix, Sol& solution, Rhs& rhs, RN const& rowNode, RTP rowTreePath, CN const& colNode, CTP colTreePath) { - Constraints::periodicBC(matrix, periodicNodes_, associations_); + std::vector> rowValues; + std::vector> colValues; + + // delete the rows (and columns) of the periodic nodes + if (symmetric_) { + for (auto const& index : periodicNodes_) { + MultiIndex row{}, col{}; + matrixWrapper(matrix.matrix()).clear_col(row, index, colValues); + matrixWrapper(matrix.matrix()).clear_row(index, col, rowValues); + } + } else { + for (auto const& index : periodicNodes_) { + MultiIndex col{}; + matrixWrapper(matrix.matrix()).clear_row(index, col, rowValues); + } + } + + // get the association or assert + auto association = [&](auto const& idx) + { + auto it = associations_.find(idx); + assert(it != associations_.end()); + return it->second; + }; + + { auto inserter = matrixInserter(matrix.matrix(), NoMatrixIndexSet{2}, false); + + // inserter {+1, -1} into left-row / right-column + for (auto const& l : periodicNodes_) { + auto&& r = association(l); + + inserter.insert(l, l, 1); + inserter.insert(l, r, -1); + if (symmetric_) { + inserter.insert(r, r, 1); + inserter.insert(r, l, -1); + } + } + + // add left-row to right-row + for (auto const& triple : rowValues) { + auto&& row = association(triple.row); + auto&& col = triple.col; + + inserter.insert(row, col, triple.value); + } + + // add left-column to right-column + if (symmetric_) { + for (auto const& triple : colValues) { + auto&& col = association(triple.col); + auto&& row = triple.row; + + inserter.insert(row, col, triple.value); + } + } + } for (auto const& a : associations_) { rhs[a.second] += rhs[a.first]; + rhs[a.first] = 0; solution[a.second] = solution[a.first]; } - Dune::Functions::interpolate(*rhs.basis(), rhs, [](auto const&) { return 0.0; }, periodicNodes_); } } // end namespace AMDiS diff --git a/src/amdis/ProblemStat.inc.hpp b/src/amdis/ProblemStat.inc.hpp index 0a942639d2b8038814a3c898a271384b12087f2d..26f31d76f5b3f090de8cdb290103989dbe52a03c 100644 --- a/src/amdis/ProblemStat.inc.hpp +++ b/src/amdis/ProblemStat.inc.hpp @@ -259,8 +259,8 @@ addDirichletBC(Predicate const& predicate, RowTreePath row, ColTreePath col, Val auto valueGridFct = makeGridFunction(values, globalBasis_->gridView()); using Range = RangeType_t; - using BcType = DirichletBC; - auto bc = std::make_unique(predicate, valueGridFct); + using BcType = DirichletBC; + auto bc = std::make_unique(predicate, valueGridFct, true); dirichletBCs_[i][j].push_back(std::move(bc)); } @@ -278,8 +278,8 @@ addDirichletBC(BoundaryType id, RowTreePath row, ColTreePath col, Values const& auto valueGridFct = makeGridFunction(values, globalBasis_->gridView()); using Range = RangeType_t; - using BcType = DirichletBC; - auto bc = std::make_unique(boundaryManager_, id, valueGridFct); + using BcType = DirichletBC; + auto bc = std::make_unique(boundaryManager_, id, valueGridFct, true); dirichletBCs_[i][j].push_back(std::move(bc)); } @@ -421,23 +421,13 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as Dune::Timer t; // 1. init matrix and rhs vector and initialize dirichlet boundary conditions - systemMatrix_->init(asmMatrix); + systemMatrix_->init(asmMatrix, dirichletBCs_, periodicBCs_); rhs_->init(asmVector); - auto localView = globalBasis_->localView(); - for_each_node(localView.tree(), [&,this](auto const& rowNode, auto rowTp) -> void { - auto rowBasis = Dune::Functions::subspaceBasis(*globalBasis_, rowTp); - for_each_node(localView.tree(), [&,this](auto const& colNode, auto colTp) -> void { - auto colBasis = Dune::Functions::subspaceBasis(*globalBasis_, colTp); - for (auto bc : dirichletBCs_[rowNode][colNode]) - bc->init(rowBasis, colBasis); - for (auto bc : periodicBCs_[rowNode][colNode]) - bc->init(rowBasis, colBasis); - }); - }); msg("{} total DOFs", globalBasis_->dimension()); // 2. traverse grid and assemble operators on the elements + auto localView = globalBasis_->localView(); for (auto const& element : elements(gridView())) { localView.bind(element); @@ -450,11 +440,11 @@ buildAfterAdapt(AdaptInfo& /*adaptInfo*/, Flag /*flag*/, bool asmMatrix, bool as } // 3. finish matrix insertion and apply dirichlet boundary conditions - systemMatrix_->finish(asmMatrix); - rhs_->finish(asmVector); + systemMatrix_->finish(); + rhs_->finish(); - for_each_node(localView.tree(), [&,this](auto const& rowNode, auto row_tp) -> void { - for_each_node(localView.tree(), [&,this](auto const& colNode, auto col_tp) -> void { + for_each_node(localView.tree(), [&](auto const& rowNode, auto row_tp) -> void { + for_each_node(localView.tree(), [&](auto const& colNode, auto col_tp) -> void { // finish boundary condition for (auto bc : dirichletBCs_[rowNode][colNode]) bc->fillBoundaryCondition(*systemMatrix_, *solution_, *rhs_, rowNode, row_tp, colNode, col_tp); diff --git a/src/amdis/ProblemStatTraits.hpp b/src/amdis/ProblemStatTraits.hpp index 6bd0c49b0cc0d38d372f1c9471103e055790fd8d..f35d0879f996d6d6b49acae2d305e402ac3c5aea 100644 --- a/src/amdis/ProblemStatTraits.hpp +++ b/src/amdis/ProblemStatTraits.hpp @@ -57,10 +57,12 @@ namespace AMDiS } // end namespace Impl /// An Exemplar for ProblemStatTraits - template + template struct DefaultProblemTraits { using GlobalBasis = GlobalBasisType; + using GridView = typename GlobalBasis::GridView; + using CoefficientType = T; }; diff --git a/src/amdis/common/CMakeLists.txt b/src/amdis/common/CMakeLists.txt index d71fa70689a09777acc3e759747d8942ff2e647d..7c5d6cbb2bdc5b92999dbe924bfb6a3e481ce68a 100644 --- a/src/amdis/common/CMakeLists.txt +++ b/src/amdis/common/CMakeLists.txt @@ -17,11 +17,10 @@ install(FILES ForEach.hpp HybridSize.hpp Index.hpp + Iterator.hpp Literals.hpp Logical.hpp Math.hpp - MultiTypeMatrix.hpp - MultiTypeVector.hpp Range.hpp Resize.hpp StaticSize.hpp diff --git a/src/amdis/common/Index.hpp b/src/amdis/common/Index.hpp index ab44576f89fe15d5fa44dc2f5e6e431cf18e75b4..dd2609a99c01e57eb5cfde19a77e7ea2493280ee 100644 --- a/src/amdis/common/Index.hpp +++ b/src/amdis/common/Index.hpp @@ -34,6 +34,13 @@ namespace AMDiS template constexpr index_t index_ = {}; + // some predefined constants + + using zero_t = index_t<0>; + static constexpr zero_t zero_ = {}; + + using one_t = index_t<1>; + static constexpr one_t one_ = {}; /// class that represents a sequence of indices template diff --git a/src/amdis/common/Iterator.hpp b/src/amdis/common/Iterator.hpp new file mode 100644 index 0000000000000000000000000000000000000000..aac66c3b4968f81e3bd8779693681c5490f195c7 --- /dev/null +++ b/src/amdis/common/Iterator.hpp @@ -0,0 +1,299 @@ +#pragma once + +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace AMDiS +{ + namespace Impl + { + // forward declaration + template + class ZipIteratorImpl; + + // forward declaration + template + struct MinimalIteratorCategory; + } + + /// The minimal iterator categories the passed iterators have in common + template + using MinimalIteratorCategory_t + = typename Impl::MinimalIteratorCategory::iterator_category...>::type; + + /// Type of a zipped iterator created by \ref zip_iterator + template + using ZipIterator + = Impl::ZipIteratorImpl, Iterators...>; + + /// Zip the passed iterator to form a new iterator that iterates in parallel through the individual elements + template + ZipIterator zip_iterator(Iterators... iterators) + { + return ZipIterator{iterators...}; + } + + + namespace Impl + { + template + struct MinimalIteratorCategory + { + using type = typename MinimalIteratorCategory::type>::type; + }; + + template + struct MinimalIteratorCategory + { + using type = std::conditional_t::value, Cat1, Cat0>; + }; + + template + struct MinimalIteratorCategory + { + using type = Cat; + }; + + + // implementation for forward iterators + template + class ZipIteratorImpl + { + public: + using iterator_category = std::forward_iterator_tag; + using value_type = std::tuple::value_type...>; + using difference_type = std::common_type_t::difference_type...>; + using reference = std::tuple::reference...>; + using pointer = void; + + template + friend class ZipIteratorImpl; + + public: + explicit ZipIteratorImpl(Iterators... iterators) + : iterators_(iterators...) + {} + + /// pre-increment all iterators + ZipIteratorImpl& operator++() + { + Tools::for_each(iterators_, [](auto& it) { ++it; }); + return *this; + } + + /// post-increment all iterators + ZipIteratorImpl operator++(int) + { + ZipIteratorImpl tmp = *this; + ++(*this); + return tmp; + } + + /// Compare all iterators by == + template + bool operator==(ZipIterator const& other) const + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + return iterators_ == other.iterators_; + } + + /// Compare all iterators by != + template + bool operator!=(ZipIterator const& other) const + { + return !(*this == other); + } + + /// Return a tuple of const references to the dereferenced iterators + auto operator*() const + { + return Tools::apply([](auto... it) { return std::tie(*it...); }, iterators_); + } + + /// Return a tuple of mutable references to the dereferenced iterators + auto operator*() + { + return Tools::apply([](auto... it) { return std::tie(*it...); }, iterators_); + } + + protected: + std::tuple iterators_; + }; + + + // Implementation for bidirectional iterators. Derived the methods from the forward iterators implementation. + template + class ZipIteratorImpl + : public ZipIteratorImpl + { + using Super = ZipIteratorImpl; + + public: + using iterator_category = std::bidirectional_iterator_tag; + using value_type = typename Super::value_type; + using difference_type = typename Super::difference_type; + using reference = typename Super::reference; + using pointer = typename Super::pointer; + + template + friend class ZipIteratorImpl; + + public: + using Super::ZipIteratorImpl; + + /// pre-decrement of all iterators + ZipIteratorImpl& operator--() + { + Tools::for_each(iterators_, [](auto& it) { --it; }); + return *this; + } + + /// post-decrement of all iterators + ZipIteratorImpl operator--(int) + { + ZipIteratorImpl tmp = *this; + --(*this); + return tmp; + } + + protected: + using Super::iterators_; + }; + + + // Implementation for random-access iterators. Derived the methods from the bidirectional iterators implementation. + template + class ZipIteratorImpl + : public ZipIteratorImpl + { + using Super = ZipIteratorImpl; + + public: + using iterator_category = std::random_access_iterator_tag; + using value_type = typename Super::value_type; + using difference_type = typename Super::difference_type; + using reference = typename Super::reference; + using pointer = typename Super::pointer; + + template + friend class ZipIteratorImpl; + + public: + using Super::ZipIteratorImpl; + + // less-than comparison of all iterators + template + bool operator<(ZipIteratorImpl const& other) const + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + return iterators_ < other.iterators_; + } + + // hreater-than comparison of all iterators + template + bool operator>(ZipIteratorImpl const& other) const + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + return iterators_ > other.iterators_; + } + + // less-than-or-equals comparison of all iterators + template + bool operator<=(ZipIteratorImpl const& other) const + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + return iterators_ <= other.iterators_; + } + + // greater-than-or-equals comparison of all iterators + template + bool operator>=(ZipIteratorImpl const& other) const + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + return iterators_ >= other.iterators_; + } + + /// Increment all iterators by shift n + ZipIteratorImpl& operator+=(difference_type n) + { + Tools::for_each(iterators_, [n](auto& it) { it += n; }); + return *this; + } + + /// Decrement all iterators by shift n + ZipIteratorImpl& operator-=(difference_type n) + { + Tools::for_each(iterators_, [n](auto& it) { it -= n; }); + return *this; + } + + /// Return shifted iterator + ZipIteratorImpl operator+(difference_type n) const + { + auto shift = [n](auto it) { return it + n; }; + return Tools::apply([shift](auto... it) { return ZipIteratorImpl(shift(it)...); }, iterators_); + } + + /// Return shifted iterator + friend ZipIteratorImpl operator+(difference_type n, ZipIteratorImpl const& zip) + { + auto shift = [n](auto it) { return n + it; }; + return Tools::apply([shift](auto... it) { return ZipIteratorImpl(shift(it)...); }, zip.iterators_); + } + + /// Return shifted iterator + ZipIteratorImpl operator-(difference_type n) const + { + auto shift = [n](auto it) { return it - n; }; + return Tools::apply([shift](auto... it) { return ZipIteratorImpl(shift(it)...); }, iterators_); + } + + + /// Return minimal difference between two zipped iterators + template + difference_type operator-(ZipIteratorImpl const& other) + { + static_assert(sizeof...(Iterators) == sizeof...(Its), ""); + auto diff = [&](auto i) -> difference_type + { + using std::get; + return get(this->iterators_) - get(other.iterators_); + }; + return Tools::apply_indices([diff](auto... i) { return Math::min(diff(i)...); }); + } + + + /// Return a tuple of const references to the dereferenced iterators + auto operator[](difference_type n) const + { + return Tools::apply([n](auto... it) { return std::tie(it[n]...); }, iterators_); + } + + /// Return a tuple of mutable references to the dereferenced iterators + auto operator[](difference_type n) + { + return Tools::apply([n](auto... it) { return std::tie(it[n]...); }, iterators_); + } + + protected: + using Super::iterators_; + }; + + } // end namespace Impl +} // end namespace AMDiS + +namespace std +{ + template + void swap(tuple lhs, tuple rhs) + { + AMDiS::Tools::for_range<0,sizeof...(Types)>([&](auto i) { swap(get(lhs), get(rhs)); }); + } +} \ No newline at end of file diff --git a/src/amdis/common/Logical.hpp b/src/amdis/common/Logical.hpp index 98e8fac4ee8f32cdc857a53475adea9f5631b366..727c493d7d611d439140bd4f7efdcfc24505d600 100644 --- a/src/amdis/common/Logical.hpp +++ b/src/amdis/common/Logical.hpp @@ -72,26 +72,16 @@ namespace AMDiS template constexpr bool_t operator!(bool_t) { return {}; } - - namespace Impl + template + struct IsEqual { - template - struct IsEqualImpl; - - template - struct IsEqualImpl - : public std::is_same, - std::integer_sequence> {}; + using type = std::is_same, std::tuple>; + }; - template - struct IsEqualImpl { enum { value = true }; }; - - } // end namespace Impl + template + using IsEqual_t = typename IsEqual::type; template - using IsEqual = Impl::IsEqualImpl; - - template - using is_one_of = or_t::value...>; + using IsEqualValue_t = IsEqual_t...>; } // end namespace AMDiS diff --git a/src/amdis/functions/Blocking.hpp b/src/amdis/functions/Blocking.hpp new file mode 100644 index 0000000000000000000000000000000000000000..859eeef120c9c61eb79837bcf83271cdfb3fd318 --- /dev/null +++ b/src/amdis/functions/Blocking.hpp @@ -0,0 +1,55 @@ +#pragma once + +#include +#include + +#include +#include + +namespace AMDiS +{ + // Lightweight representation of (hierarchic) block structure + namespace tag + { + /// \brief Blocking structure corresponding to IndexMergingStrategies + /// \ref FlatLexicographic and \ref FlatInterleaved + struct flat + { + template + flat operator[](const Index&) const { return {}; } + }; + + /// \brief Blocking structure corresponding to IndexMergingStrategy + /// \ref LeafBlockedInterleaved + template + struct leaf_blocked + { + template + flat operator[](const Index&) const { return {}; } + }; + + /// \brief Blocking structure corresponding to IndexMergingStrategy + /// \ref BlockedLexicographic + template + struct blocked + { + template + auto operator[](index_t) const { return std::get(t); } + Tag0 operator[](std::size_t /*i*/) const { return {}; } + + std::tuple t; + }; + + } // end namespace tag + + // forward declaration + template + struct Blocking; + + /// \brief Extract blocking structure of GlobalBasis + template + using Blocking_t = typename Blocking::type; + +} // end namespace AMDiS + +#include diff --git a/src/amdis/functions/Blocking.inc.hpp b/src/amdis/functions/Blocking.inc.hpp new file mode 100644 index 0000000000000000000000000000000000000000..297399ecda4f59a3a7c0aa2162d95418778ad9c2 --- /dev/null +++ b/src/amdis/functions/Blocking.inc.hpp @@ -0,0 +1,119 @@ +#pragma once + +#include +#include + +#include +#include +#include +#include + +#include +#include +#include + +namespace AMDiS +{ + namespace Impl + { + // Handle non-leaf nodes + template