diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000000000000000000000000000000000000..442ac76345bc8676ac9ff97ebddd226b829c9542
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "externals/fmt"]
+ path = externals/fmt
+ url = https://github.com/fmtlib/fmt.git
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 1cb8e11196d5fe19bd37d25d880c4a789b00c143..96780b564c01ae7e5eaf11d4036a83cc8ad94ee1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -19,6 +19,7 @@ add_subdirectory("test")
add_subdirectory("examples" EXCLUDE_FROM_ALL)
add_subdirectory("doc")
add_subdirectory("cmake/modules")
+add_subdirectory("externals" EXCLUDE_FROM_ALL)
target_include_directories(amdis PUBLIC
$<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/src>)
diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt
index ab96c269b68d746d67d4b4eda919c38d785fd312..d799f6d14ed1e55fc0079930b91b678d0c568e21 100644
--- a/examples/CMakeLists.txt
+++ b/examples/CMakeLists.txt
@@ -6,7 +6,7 @@ foreach(project ${projects2d})
add_executable(${project}.2d ${project}.cc)
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project}.2d)
target_link_dune_default_libraries(${project}.2d)
- target_link_libraries(${project}.2d "amdis")
+ target_link_libraries(${project}.2d amdis fmt)
target_compile_definitions(${project}.2d PRIVATE AMDIS_DIM=2 AMDIS_DOW=2)
add_dependencies(examples ${project}.2d)
endforeach()
@@ -17,7 +17,7 @@ foreach(project ${projects3d})
add_executable(${project}.3d ${project}.cc)
add_dune_alberta_flags(GRIDDIM 3 WORLDDIM 3 ${project}.3d)
target_link_dune_default_libraries(${project}.3d)
- target_link_libraries(${project}.3d "amdis")
+ target_link_libraries(${project}.3d amdis fmt)
target_compile_definitions(${project}.3d PRIVATE AMDIS_DIM=3 AMDIS_DOW=3)
add_dependencies(examples ${project}.3d)
endforeach()
\ No newline at end of file
diff --git a/examples/ellipt.cc b/examples/ellipt.cc
index 6baf47bfb8cbba1adf2a4ffc10a5cb46f742a29f..0414e633a00c14d0c0ee320d7728b90583257d3e 100644
--- a/examples/ellipt.cc
+++ b/examples/ellipt.cc
@@ -3,55 +3,95 @@
#include <iostream>
+#include <fmt/core.h>
+
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/Operators.hpp>
#include <amdis/common/Literals.hpp>
+#include <amdis/gridfunctions/Integrate.hpp>
using namespace AMDiS;
+using namespace Dune::Indices;
// 1 component with polynomial degree 1
-//using Grid = Dune::AlbertaGrid<AMDIS_DIM, AMDIS_DOW>;
-using ElliptParam = YaspGridBasis<AMDIS_DIM, 1>;
-using ElliptProblem = ProblemStat<ElliptParam>;
+using Param = YaspGridBasis<AMDIS_DIM, 1>;
+using ElliptProblem = ProblemStat<Param>;
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
- using namespace Dune::Indices;
+ int numLevels = AMDIS_DIM == 2 ? 8 : 5;
+ if (argc > 2)
+ numLevels = std::atoi(argv[2]);
+
+ auto f = [](auto const& x) {
+ double r2 = dot(x,x);
+ double ux = std::exp(-10.0 * r2);
+ return -(400.0 * r2 - 20.0 * x.size()) * ux;
+ };
+ auto g = [](auto const& x){ return std::exp(-10.0 * dot(x,x)); };
+ auto grad_g = [](auto const& x) -> FieldMatrix<double,1,AMDIS_DIM> {
+ return {-20.0 * std::exp(-10.0 * dot(x,x)) * x};
+ };
ElliptProblem prob("ellipt");
prob.initialize(INIT_ALL);
- AdaptInfo adaptInfo("adapt");
-
auto opL = makeOperator(tag::gradtest_gradtrial{}, 1.0);
prob.addMatrixOperator(opL, _0, _0);
- auto opForce = makeOperator(tag::test{}, [](auto const& x) { return -1.0; }, 0);
+ auto opForce = makeOperator(tag::test{}, f, 6);
prob.addVectorOperator(opForce, _0);
- auto opForce2 = makeOperator(tag::test{}, [](auto const& x) { return -2.0; }, 0);
- prob.addVectorOperator(BoundaryType{0}, opForce2, _0);
-
// set boundary condition
- auto predicate = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8; }; // define boundary
- auto dbcValues = [](auto const& x){ return 0.0; }; // set value
- prob.addDirichletBC(predicate, _0, _0, dbcValues);
+ auto boundary = [](auto const& x){ return x[0] < 1.e-8 || x[1] < 1.e-8 || x[0] > 1.0-1.e-8 || x[1] > 1.0-1.e-8; };
+ prob.addDirichletBC(boundary, _0, _0, g);
+
+ AdaptInfo adaptInfo("adapt");
- prob.buildAfterCoarsen(adaptInfo, Flag(0));
+ std::vector<double> errL2; errL2.reserve(numLevels);
+ std::vector<double> errH1; errH1.reserve(numLevels);
+ std::vector<double> widths; widths.reserve(numLevels);
+ for (int i = 0; i < numLevels; ++i) {
+ prob.getGrid()->globalRefine(1);
+ auto gridView = prob.getGlobalBasis()->gridView();
+
+ double h = 0;
+ for (auto const& e : edges(gridView))
+ h = std::max(h, e.geometry().volume());
+ widths.push_back(h);
+
+ prob.buildAfterCoarsen(adaptInfo, Flag(0));
+ prob.solve(adaptInfo);
+
+ double errorL2 = integrate(sqr(g - prob.getSolution(_0)), gridView, 6);
+ errL2.push_back(std::sqrt(errorL2));
+ double errorH1 = errorL2 + integrate(unary_dot(grad_g - gradientAtQP(prob.getSolution(_0))), gridView, 6);
+ errH1.push_back(std::sqrt(errorH1));
+
+#if WRITE_FILES
+ Dune::VTKWriter<typename ElliptProblem::GridView> vtkWriter(prob.getGlobalBasis()->gridView());
+ vtkWriter.addVertexData(prob.getSolution(_0), Dune::VTK::FieldInfo("u", Dune::VTK::FieldInfo::Type::scalar, 1));
+ vtkWriter.write("u_" + std::to_string(i));
+#endif
+ }
- // write matrix to file
- if (Parameters::get<int>("elliptMesh->global refinements").value_or(0) < 5) {
- mtl::io::matrix_market_ostream out("matrix.mtx");
- out << prob.getSystemMatrix()->getMatrix();
+ std::cout << "\n";
+ fmt::print("{:5} | {:12} | {:12} | {:12} | {:12} | {:12}\n",
+ "level", "h", "|u - u*|_L2","|u - u*|_H1","eoc_L2","eoc_H1");
+ fmt::print("{0:->7}{0:->15}{0:->15}{0:->15}{0:->15}{1:->14}","+","\n");
- std::cout << prob.getSystemMatrix()->getMatrix() << '\n';
+ std::vector<double> eocL2(numLevels, 0.0), eocH1(numLevels, 0.0);
+ for (int i = 1; i < numLevels; ++i) {
+ eocL2[i] = std::log(errL2[i]/errL2[i-1]) / std::log(widths[i]/widths[i-1]);
+ eocH1[i] = std::log(errH1[i]/errH1[i-1]) / std::log(widths[i]/widths[i-1]);
}
- prob.solve(adaptInfo);
- prob.writeFiles(adaptInfo, true);
+ for (int i = 0; i < numLevels; ++i)
+ fmt::print("{:<5} | {:<12} | {:<12} | {:<12} | {:<12} | {:<12}\n",
+ i+1, widths[i], errL2[i], errH1[i], eocL2[i], eocH1[i]);
AMDiS::finalize();
return 0;
diff --git a/examples/init/ellipt.dat.2d b/examples/init/ellipt.dat.2d
index 37ba6ea98163b277213cc64c1a5afe3c1f823610..f12565729e95e06f7902cf81744b29450c3445e1 100644
--- a/examples/init/ellipt.dat.2d
+++ b/examples/init/ellipt.dat.2d
@@ -1,13 +1,10 @@
-dimension of world: 2
-
-elliptMesh->macro file name: ./macro/macro.stand.2d
-elliptMesh->global refinements: 5
+elliptMesh->global refinements: 0
ellipt->mesh: elliptMesh
-ellipt->solver->name: cg
+ellipt->solver->name: bicgstab
ellipt->solver->max iteration: 1000
-ellipt->solver->tolerance: 1.e-8
+ellipt->solver->relative tolerance: 1.e-8
ellipt->solver->info: 10
ellipt->solver->left precon: diag
ellipt->solver->right precon: no
diff --git a/examples/init/ellipt.dat.3d b/examples/init/ellipt.dat.3d
new file mode 100644
index 0000000000000000000000000000000000000000..c05c023447ac2718894383a4d5df234cfbe4dade
--- /dev/null
+++ b/examples/init/ellipt.dat.3d
@@ -0,0 +1,14 @@
+elliptMesh->global refinements: 0
+
+ellipt->mesh: elliptMesh
+
+ellipt->solver->name: bicgstab
+ellipt->solver->max iteration: 1000
+ellipt->solver->relative tolerance: 1.e-8
+ellipt->solver->info: 10
+ellipt->solver->left precon: diag
+ellipt->solver->right precon: no
+
+ellipt->output[0]->filename: ellipt.3d
+ellipt->output[0]->name: u
+ellipt->output[0]->output directory: ./output
diff --git a/externals/CMakeLists.txt b/externals/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..763f19e11d62943509a44746ba6c70ad04bd390e
--- /dev/null
+++ b/externals/CMakeLists.txt
@@ -0,0 +1 @@
+add_subdirectory(fmt)
\ No newline at end of file
diff --git a/externals/fmt b/externals/fmt
new file mode 160000
index 0000000000000000000000000000000000000000..c6d9730ddbfaa5dbefd5d2b48ef047ad362157f4
--- /dev/null
+++ b/externals/fmt
@@ -0,0 +1 @@
+Subproject commit c6d9730ddbfaa5dbefd5d2b48ef047ad362157f4
diff --git a/src/amdis/DirichletBC.hpp b/src/amdis/DirichletBC.hpp
index 1729ec65a3f51a56f52450b108d924f236730ad5..d456189e62747eff3bad6265f70af886c41c4e73 100644
--- a/src/amdis/DirichletBC.hpp
+++ b/src/amdis/DirichletBC.hpp
@@ -80,6 +80,7 @@ namespace AMDiS
// subtract columns of dirichlet nodes from rhs
// for (auto const& triplet : columns)
// rhs[triplet.row] -= triplet.value * solution[triplet.col];
+ initialized_ = false;
}
protected:
diff --git a/src/amdis/GridFunctionOperator.hpp b/src/amdis/GridFunctionOperator.hpp
index 70a12642d81e1fd311b8957c6833bf0e751378f2..318506fa12473ce9c561258b8a9f781d4cbb2975 100644
--- a/src/amdis/GridFunctionOperator.hpp
+++ b/src/amdis/GridFunctionOperator.hpp
@@ -61,24 +61,9 @@ namespace AMDiS
**/
GridFunctionOperatorBase(GridFunction const& gridFct, int termOrder)
: gridFct_(gridFct)
- , localFct_(localFunction(gridFct_))
, termOrder_(termOrder)
{}
- // Copy constructor
- GridFunctionOperatorBase(GridFunctionOperatorBase const& that)
- : gridFct_(that.gridFct_)
- , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
- , termOrder_(that.termOrder_)
- {}
-
- // Move constructor
- GridFunctionOperatorBase(GridFunctionOperatorBase&& that)
- : gridFct_(std::move(that.gridFct_))
- , localFct_(localFunction(gridFct_)) // recreate localFct_ on new gridFct_
- , termOrder_(std::move(that.termOrder_))
- {}
-
/// \brief Binds operator to `element` and `geometry`.
/**
* Binding an operator to the currently visited element in grid traversal.
@@ -91,21 +76,24 @@ namespace AMDiS
void bind_impl(Element const& element, Geometry const& geometry)
{
assert( bool(quadFactory_) );
- localFct_.bind(element);
- quadFactory_->bind(localFct_);
+ localFct_.emplace(localFunction(gridFct_));
+ localFct_->bind(element);
+ quadFactory_->bind(localFct_.value());
}
/// Unbinds operator from element.
void unbind_impl()
{
- localFct_.unbind();
+ localFct_->unbind();
+ localFct_.reset();
}
/// Create a quadrature factory from a PreQuadratureFactory, e.g. class derived from \ref QuadratureFactory
template <class PreQuadFactory>
void setQuadFactory(PreQuadFactory const& pre)
{
- quadFactory_ = makeQuadratureFactoryPtr<typename Geometry::ctype, LocalContext::mydimension, LocalFunction>(pre);
+ using ctype = typename Geometry::ctype;
+ quadFactory_.emplace(makeQuadratureFactory<ctype, LocalContext::mydimension, LocalFunction>(pre));
}
protected:
@@ -113,7 +101,7 @@ namespace AMDiS
auto coefficient(LocalCoordinate const& local) const
{
assert( this->bound_ );
- return localFct_(local);
+ return (*localFct_)(local);
}
/// Create a quadrature rule using the \ref QuadratureFactory by calculating the
@@ -131,10 +119,10 @@ namespace AMDiS
GridFunction gridFct_;
/// localFunction associated with gridFunction. Mus be updated whenever gridFunction changes.
- LocalFunction localFct_;
+ Dune::Std::optional<LocalFunction> localFct_;
/// Assign each element type a quadrature rule
- std::shared_ptr<QuadFactory> quadFactory_ = std::make_shared<QuadFactory>();
+ Dune::Std::optional<QuadFactory> quadFactory_;
/// the derivative order of this operator (in {0, 1, 2})
int termOrder_ = 0;
diff --git a/src/amdis/common/Math.hpp b/src/amdis/common/Math.hpp
index f3a7f3b3ae6718fede5558194ab59cc22e422559..555d1585c0cbb0cb97bd83e727f7afa681845834 100644
--- a/src/amdis/common/Math.hpp
+++ b/src/amdis/common/Math.hpp
@@ -118,7 +118,7 @@ namespace AMDiS
}
template <class T>
- constexpr T threshold = T(1.e-16); //Math::sqr(std::numeric_limits<T>::epsilon());
+ constexpr T threshold = T(1.e-18L); //Math::sqr(std::numeric_limits<T>::epsilon());
/// Calculates factorial of i
diff --git a/src/amdis/linear_algebra/LinearSolverInterface.hpp b/src/amdis/linear_algebra/LinearSolverInterface.hpp
index f308e1d2c5200cae184a45a85348cfb2e00c7254..5d1534f9f522db80170bdd3ee36bea85d25d56cd 100644
--- a/src/amdis/linear_algebra/LinearSolverInterface.hpp
+++ b/src/amdis/linear_algebra/LinearSolverInterface.hpp
@@ -45,6 +45,7 @@ namespace AMDiS
void solve(Matrix const& A, VectorX& x, VectorB const& b,
SolverInfo& solverInfo)
{
+ assert( num_cols(A) == size(x) && num_rows(A) == size(b) );
solveImpl(A, x, b, solverInfo);
}
diff --git a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
index efadcc319d8618f27bb854331f083a08bbbe4730..dba705f9dba0809107058fb1ea8e3e12dd71ccb3 100644
--- a/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFMatrix.hpp
@@ -146,7 +146,9 @@ namespace AMDiS
test_exit(!insertionMode && !inserter, "Matrix already in insertion mode!");
calculateSlotSize();
- matrix->change_dim(rowFeSpace.size(), colFeSpace.size());
+ matrix->change_dim(rowFeSpace.dimension(), colFeSpace.dimension());
+ test_exit(num_rows(*matrix) == rowFeSpace.dimension() && num_cols(*matrix) == colFeSpace.dimension(),
+ "Wrong dimensions in matrix");
if (prepareForInsertion) {
set_to_zero(*matrix);
inserter = new Inserter(*matrix, slot_size);
@@ -281,7 +283,7 @@ namespace AMDiS
private:
/// The minimal number of nonzeros per row
- static constexpr int MIN_NNZ_PER_ROW = 20;
+ static constexpr int MIN_NNZ_PER_ROW = 50;
/// The finite element space / basis of the row
RowFeSpace const& rowFeSpace;
@@ -302,7 +304,7 @@ namespace AMDiS
bool insertionMode = false;
/// The size of the slots used to insert values per row
- int slot_size = 20;
+ int slot_size = MIN_NNZ_PER_ROW;
// friend class declarations
template <class, class> friend class SystemMatrix;
diff --git a/src/amdis/linear_algebra/mtl/DOFVector.hpp b/src/amdis/linear_algebra/mtl/DOFVector.hpp
index b66d6ab5cf2c7bb71f46f6c601d24dc1e1c56e3c..6cfb293aac137b7c371a73642be289c0f64e54cb 100644
--- a/src/amdis/linear_algebra/mtl/DOFVector.hpp
+++ b/src/amdis/linear_algebra/mtl/DOFVector.hpp
@@ -45,7 +45,6 @@ namespace AMDiS
, vector(ClonablePtr<BaseVector>::make())
{
compress();
- *vector = ValueType(0);
}
/// Constructor. Takes reference to existing BaseVector
@@ -89,8 +88,10 @@ namespace AMDiS
/// Resize the \ref vector to the size of the \ref feSpace.
void compress()
{
- if (num_rows(*vector) != getSize())
+ if (num_rows(*vector) != getSize()) {
resize(getSize());
+ *vector = ValueType(0);
+ }
}
template <class SizeInfo>