diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index fa2600207858d788526f695d0d09861ebc4485fb..15799aaa6707c156c45780e88c19aee2cea29826 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -31,23 +31,4 @@ 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 () - -add_custom_target(ellipt) -foreach(_D 2 3) - add_amdis_executable(NAME ellipt.yaspgrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D}) - target_compile_definitions(ellipt.yaspgrid.${_D}d PRIVATE "GRID_ID=0") - add_dependencies(ellipt ellipt.yaspgrid.${_D}d) - - add_amdis_executable(NAME ellipt.albertagrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D} ALBERTA_GRID) - target_compile_definitions(ellipt.albertagrid.${_D}d PRIVATE "GRID_ID=1") - add_dependencies(ellipt ellipt.albertagrid.${_D}d) - - add_amdis_executable(NAME ellipt.uggrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D}) - target_compile_definitions(ellipt.uggrid.${_D}d PRIVATE "GRID_ID=2") - add_dependencies(ellipt ellipt.uggrid.${_D}d) - - add_amdis_executable(NAME ellipt.alugrid.${_D}d SOURCES ellipt.cc DIM ${_D} DOW ${_D}) - target_compile_definitions(ellipt.alugrid.${_D}d PRIVATE "GRID_ID=3") - add_dependencies(ellipt ellipt.alugrid.${_D}d) -endforeach() +endif () \ No newline at end of file diff --git a/examples/ellipt.cc b/examples/ellipt.cc index c96a319ec19a103c0b6951e0d9a51923574ab100..411ef19500b5343088b9b71053dd7c6ad66b6784 100644 --- a/examples/ellipt.cc +++ b/examples/ellipt.cc @@ -1,32 +1,26 @@ +#ifdef HAVE_CONFIG_H +#include "config.h" +#endif +#include <iostream> + #include <amdis/AMDiS.hpp> #include <amdis/Integrate.hpp> #include <amdis/LocalOperators.hpp> #include <amdis/ProblemStat.hpp> #include <amdis/common/Literals.hpp> -#if GRID_ID == 1 -#include <dune/grid/albertagrid.hh> -using Grid = Dune::AlbertaGrid<GRIDDIM,WORLDDIM>; -#elif GRID_ID == 2 -#include <dune/grid/uggrid.hh> -using Grid = Dune::UGGrid<GRIDDIM>; -#elif GRID_ID == 3 -#include <dune/alugrid/grid.hh> -using Grid = Dune::ALUGrid<GRIDDIM,WORLDDIM,Dune::simplex,Dune::conforming>; -#else -#include <dune/grid/yaspgrid.hh> -using Grid = Dune::YaspGrid<GRIDDIM>; -#endif - using namespace AMDiS; using namespace Dune::Indices; +// 1 component with polynomial degree 1 +using Param = YaspGridBasis<GRIDDIM, 2>; +using ElliptProblem = ProblemStat<Param>; + int main(int argc, char** argv) { Environment env(argc, argv); - Dune::Timer t; - int numLevels = GRIDDIM == 2 ? 6 : 4; + int numLevels = GRIDDIM == 2 ? 8 : 5; if (argc > 2) numLevels = std::atoi(argv[2]); @@ -40,8 +34,7 @@ int main(int argc, char** argv) return {-20.0 * std::exp(-10.0 * dot(x,x)) * x}; }; - using Param = LagrangeBasis<Grid, 2>; - ProblemStat<Param> prob("ellipt"); + ElliptProblem prob("ellipt"); prob.initialize(INIT_ALL); auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0); @@ -51,49 +44,40 @@ int main(int argc, char** argv) prob.addVectorOperator(opForce, 0); // set boundary condition -#if GRIDDIM == 2 - prob.boundaryManager()->setBoxBoundary({1,1,1,1}); -#elif GRIDDIM == 3 - prob.boundaryManager()->setBoxBoundary({1,1,1,1,1,1}); -#endif - prob.addDirichletBC(1, 0, 0, g); + prob.addDirichletBC(BoundaryType{1}, 0, 0, g); AdaptInfo adaptInfo("adapt"); std::vector<double> errL2; errL2.reserve(numLevels); std::vector<double> errH1; errH1.reserve(numLevels); std::vector<double> widths; widths.reserve(numLevels); - for (int l = 0; l < numLevels; ++l) { - prob.globalRefine(1); + for (int i = 0; i < numLevels; ++i) { + prob.grid()->globalRefine(1); + auto gridView = prob.gridView(); double h = 0; - for (auto const& e : elements(prob.gridView(), Dune::Partitions::interior)) { - auto refElem = Dune::referenceElement<double,GRIDDIM>(e.type()); - auto geo = e.geometry(); - for (int i = 0; i < refElem.size(GRIDDIM-1); ++i) { // edges - auto v0 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,0,GRIDDIM), GRIDDIM)); - auto v1 = geo.global(refElem.position(refElem.subEntity(i,GRIDDIM-1,1,GRIDDIM), GRIDDIM)); - h = std::max(h, (v0 - v1).two_norm()); - } - } + for (auto const& e : edges(gridView)) + h = std::max(h, e.geometry().volume()); widths.push_back(h); + prob.globalBasis()->update(gridView); + prob.solutionVector()->compress(); prob.assemble(adaptInfo); prob.solve(adaptInfo); - double errorL2 = integrate(sqr(g - prob.solution(0)), prob.gridView(), 6); + double errorL2 = integrate(sqr(g - prob.solution(0)), gridView, 6); errL2.push_back(std::sqrt(errorL2)); - double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), prob.gridView(), 6); + double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.solution(0))), gridView, 6); errH1.push_back(std::sqrt(errorH1)); #if WRITE_FILES Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(gridView); vtkWriter.addVertexData(prob.solution(0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1)); - vtkWriter.write("u_" + std::to_string(l)); + vtkWriter.write("u_" + std::to_string(i)); #endif } - msg(""); + std::cout << "\n"; msg("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}", "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1"); msg_("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n"); @@ -108,6 +92,5 @@ int main(int argc, char** argv) msg("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}", i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]); - msg("elapsed time: {} seconds", t.elapsed()); return 0; } diff --git a/examples/init/boundary.dat.2d b/examples/init/boundary.dat.2d index ab66e76c319bb7ae16203428ff6a2f4595d8f815..23f56aeb32aec0b735cbd525033d66806802add7 100644 --- a/examples/init/boundary.dat.2d +++ b/examples/init/boundary.dat.2d @@ -1,6 +1,6 @@ -elliptMesh->global refinements: 3 +mesh->global refinements: 3 -ellipt->mesh: elliptMesh +ellipt->mesh: mesh ellipt->solver: bicgstab ellipt->solver->max iteration: 10000 diff --git a/examples/init/cahn_hilliard.dat.2d b/examples/init/cahn_hilliard.dat.2d index 71e5527a59287547e7d8e8e197f20b532fa4bd0c..56ef77a4a518355550ca9c308f428d8caa94f335 100644 --- a/examples/init/cahn_hilliard.dat.2d +++ b/examples/init/cahn_hilliard.dat.2d @@ -14,7 +14,7 @@ refinement->bulk: 2 ch->mesh: chMesh ch->solver: direct ch->solver->max iteration: 1000 -ch->solver->absolute tolerance: 1e-6 +ch->solver->relative tolerance: 1e-6 ch->solver->break if tolerance not reached: 1 ch->solver->info: 2 diff --git a/examples/init/ellipt.dat.2d b/examples/init/ellipt.dat.2d index 68fd6e3c5c9dc1e6b594a3e91b8a424399b7863c..bc4dcc23e22a85e85ac2fd9265235932896e39c4 100644 --- a/examples/init/ellipt.dat.2d +++ b/examples/init/ellipt.dat.2d @@ -1,7 +1,6 @@ elliptMesh->global refinements: 0 elliptMesh->num cells: 8 8 elliptMesh->overlap: 1 -% elliptMesh->structured: simplex ellipt->mesh: elliptMesh ellipt->symmetry: spd diff --git a/examples/init/ellipt.dat.3d b/examples/init/ellipt.dat.3d index 90a3e25b58100565f0cbe86575707cf76241cf29..ad5cbc12a751622277cb1b5e654eda8b4f8b2141 100644 --- a/examples/init/ellipt.dat.3d +++ b/examples/init/ellipt.dat.3d @@ -1,7 +1,6 @@ elliptMesh->global refinements: 0 elliptMesh->num cells: 8 8 8 elliptMesh->overlap: 1 -% elliptMesh->structured: simplex ellipt->mesh: elliptMesh ellipt->symmetry: spd diff --git a/examples/init/heat.dat.2d b/examples/init/heat.dat.2d index b58e13e031b6221fadafc57f5c7ad35f6578bbe9..4db5cd5e3063245bf794d56eb917e2bc01af6d4c 100644 --- a/examples/init/heat.dat.2d +++ b/examples/init/heat.dat.2d @@ -5,6 +5,7 @@ heatMesh->num cells: 8 8 heat->mesh: heatMesh heat->names: u + heat->solver: pcg heat->solver->max iteration: 1000 heat->solver->absolute tolerance: 1e-6 diff --git a/src/amdis/BiLinearForm.hpp b/src/amdis/BiLinearForm.hpp index 992afd5f3fe75c027fff81d59ff15b24338e9004..794a8ea2c8586544b70b2cdc8c0d432e963c7de4 100644 --- a/src/amdis/BiLinearForm.hpp +++ b/src/amdis/BiLinearForm.hpp @@ -24,7 +24,8 @@ namespace AMDiS : public MatrixBase<RB,CB,MatrixBackend<BackendTraits<RB,T>>> { using Self = BiLinearForm; - using Backend = MatrixBackend<BackendTraits<RB,T>>; + using Traits = BackendTraits<RB,T>; + using Backend = MatrixBackend<Traits>; using Super = MatrixBase<RB,CB,Backend>; public: diff --git a/src/amdis/linearalgebra/MatrixBase.hpp b/src/amdis/linearalgebra/MatrixBase.hpp index bf5d5c5ddf6a857e3448c79cd43d833df2c93acc..f0edfcd0aee77931dde65212133d6d5721d63b75 100644 --- a/src/amdis/linearalgebra/MatrixBase.hpp +++ b/src/amdis/linearalgebra/MatrixBase.hpp @@ -8,6 +8,7 @@ #include <amdis/common/Concepts.hpp> #include <amdis/common/TypeTraits.hpp> #include <amdis/functions/NodeIndices.hpp> +#include <amdis/linearalgebra/SymmetryStructure.hpp> #include <amdis/typetree/MultiIndex.hpp> namespace AMDiS diff --git a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp index 2b81e07c43d1f37ef33e67e2ed2457bd9a96ebf8..459d5203c2f408b2a0d87789f7e580b6f415d415 100644 --- a/src/amdis/linearalgebra/eigen/MatrixBackend.hpp +++ b/src/amdis/linearalgebra/eigen/MatrixBackend.hpp @@ -8,6 +8,7 @@ #include <dune/common/timer.hh> #include <amdis/Output.hpp> +#include <amdis/linearalgebra/SymmetryStructure.hpp> namespace AMDiS { diff --git a/src/amdis/linearalgebra/istl/MatrixBackend.hpp b/src/amdis/linearalgebra/istl/MatrixBackend.hpp index 7da410e61566752434599f750248435f0ca7241c..4d0b4a65545871e7197bbed4ffb203658322e969 100644 --- a/src/amdis/linearalgebra/istl/MatrixBackend.hpp +++ b/src/amdis/linearalgebra/istl/MatrixBackend.hpp @@ -7,6 +7,7 @@ #include <dune/istl/matrixindexset.hh> #include <amdis/Output.hpp> +#include <amdis/linearalgebra/SymmetryStructure.hpp> namespace AMDiS { diff --git a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp index d1f5e19a8ff6dbaad3f25f858f43aea930df378a..75820f5412a70fa5a7e0b8bfd64319414ef02def 100644 --- a/src/amdis/linearalgebra/mtl/MatrixBackend.hpp +++ b/src/amdis/linearalgebra/mtl/MatrixBackend.hpp @@ -10,6 +10,7 @@ #include <boost/numeric/mtl/utility/range_wrapper.hpp> #include <amdis/Output.hpp> +#include <amdis/linearalgebra/SymmetryStructure.hpp> namespace AMDiS {