Commit 493a7246 authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

Added fmt as output formatting library

parent c57c5580
......@@ -66,7 +66,7 @@ namespace AMDiS
value_type const& operator[](size_type i) const
{
test_exit_dbg(i < vector_.size(),
"Index ", i, " out of range [0,", vector_.size(), ")" );
"Index {} out of range [0,{})", i, vector_.size());
return vector_[i];
}
......@@ -74,7 +74,7 @@ namespace AMDiS
value_type& operator[](size_type i)
{
test_exit_dbg(i < vector_.size(),
"Index ", i, " out of range [0,", vector_.size(), ")" );
"Index {} out of range [0,{})", i, vector_.size());
return vector_[i];
}
......
......@@ -61,7 +61,7 @@ namespace AMDiS
{
test_exit_dbg(inserter_, "Inserter not initilized!");
test_exit_dbg(r < num_rows(matrix_) && c < num_cols(matrix_),
"Indices out of range [0,", num_rows(matrix_), ")x[0,", num_cols(matrix_), ")" );
"Indices out of range [0,{})x[0,{})", num_rows(matrix_), num_cols(matrix_));
(*inserter_)[r][c] += value;
}
......
......@@ -54,7 +54,7 @@ namespace AMDiS
value_type const& operator[](size_type i) const
{
test_exit_dbg(i < mtl::vec::size(vector_),
"Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
"Index {} out of range [0,{})", i, mtl::vec::size(vector_));
return vector_[i];
}
......@@ -62,7 +62,7 @@ namespace AMDiS
value_type& operator[](size_type i)
{
test_exit_dbg(i < mtl::vec::size(vector_),
"Index ", i, " out of range [0,", mtl::vec::size(vector_), ")" );
"Index {} out of range [0,{})", i, mtl::vec::size(vector_));
return vector_[i];
}
......
......@@ -60,7 +60,7 @@ namespace AMDiS
try {
code = (*solver_)(x, b);
} catch (mtl::mat::umfpack::error& e) {
error_exit("UMFPACK_ERROR(solve, ", e.code, ") = ", e.what());
error_exit("UMFPACK_ERROR(solve, {}) = {}", e.code, e.what());
}
auto r = Vector(b);
......@@ -106,7 +106,7 @@ namespace AMDiS
try {
Super::solver_.reset(new SolverType(fullMatrix, Super::symmetricStrategy_, Super::allocInit_));
} catch (mtl::mat::umfpack::error const& e) {
error_exit("UMFPACK_ERROR(factorize, ", e.code, ") = ", e.what());
error_exit("UMFPACK_ERROR(factorize, {}) = {}", e.code, e.what());
}
}
};
......
# maybe unncessary (?)
dune_add_test(NAME test_duneamdis
TARGET duneamdis
COMPILE_ONLY)
# find_package(GTest REQUIRED)
set(projects "test1" "test2")
foreach(project ${projects})
dune_add_test(NAME ${project}
SOURCES ${project}.cc
LINK_LIBRARIES duneamdis ${GTEST_BOTH_LIBRARIES}
COMPILE_DEFINITIONS "DIM=2;DOW=2"
CMD_ARGS "init/test.json.2d")
add_dune_alberta_flags(GRIDDIM 2 WORLDDIM 2 ${project})
target_include_directories(${project} PRIVATE ${GTEST_INCLUDE_DIRS})
set_tests_properties(${project} PROPERTIES WORKING_DIRECTORY ${CMAKE_SOURCE_DIR})
endforeach()
\ No newline at end of file
#pragma once
#include <nanoflann.hpp>
#include <cstdlib>
#include <iostream>
namespace AMDiS
{
template <class FeSpace>
struct KDTreeMeshAdaptor
{
using MeshView = typename FeSpace::GridView;
using SizeType = typename FeSpace::LocalIndexSet::size_type;
using ValueType = double;
static constexpr dimension = MeshView::dimension;
static constexpr dimensionworld = MeshView::dimensionworld;
using Codim0 = typename GridView::template Codim<0>;
using PointType = typename Codim0::Geometry::GlobalCoordinate;
using VectorOfPointsType = std::vector<PointType> ;
using DataType = typename MeshView::Grid::EntitySeed;
using VectorOfDataType = std::vector<DataType>;
using Distance = nanoflann::metric_L2_Simple;
using Self = KDTreeMeshAdaptor<FeSpace>;
using Metric = typename Distance::template traits<ValueType, Self>::distance_t;
using IndexType = nanoflann::KDTreeSingleIndexAdaptor<Metric, Self, dimensionworld, SizeType>;
public:
unique_ptr<IndexType> index; //! The kd-tree index for the user to call its methods as usual with any other FLANN index.
FeSpace const& feSpace;
VectorOfPointsType_ m_points;
VectorOfDataType_ m_data;
std::map<std::hash<DataType>, DataType> shapeValues;
public:
/// Constructor: takes a const ref to the vector of vectors object with the data points
KDTreeMeshAdaptor(FeSpace const& feSpace_,
const int leaf_max_size = 10)
: feSpace(feSpace_)
, localView(feSpace.localView())
{
init();
index = make_unique<IndexType>(dimensionworld, *this /* adaptor */,
nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
index->buildIndex();
}
/// Destructor.
~KDTreeMeshAdaptor() {}
void init()
{
for (auto const& element : elements(feSpace.gridView())) {
m_points.push_back(element.geometry().center());
m_data.push_back(element.seed());
}
}
void reinit(const int leaf_max_size = 10)
{
m_points.clear();
m_data.clear();
init();
index.reset(new IndexType(dimensionworld, *this /* adaptor */,
nanoflann::KDTreeSingleIndexAdaptorParams(leaf_max_size, dimensionworld));
index->buildIndex();
}
/** Query for the \a num_closest closest points to a given point (entered as query_point[0:dim-1]).
* Note that this is a short-cut method for index->findNeighbors().
* The user can also call index->... methods as desired.
* \note nChecks_IGNORED is ignored but kept for compatibility with the original FLANN interface.
*/
inline void query(const ValueType *query_point,
const std::size_t num_closest,
SizeType *out_indices,
ValueType *out_distances_sq,
const int nChecks_IGNORED = 10) const
{
nanoflann::KNNResultSet<ValueType, SizeType> resultSet(num_closest);
resultSet.init(out_indices, out_distances_sq);
index->findNeighbors(resultSet, query_point, nanoflann::SearchParams());
}
/** @name Interface expected by KDTreeSingleIndexAdaptor
* @{ */
Self const& derived() const
{
return *this;
}
Self& derived()
{
return *this;
}
// Must return the number of data points
inline std::size_t kdtree_get_point_count() const
{
return m_points.size();
}
// Returns the distance between the vector "p1[0:size-1]" and the data point with index "idx_p2" stored in the class:
inline ValueType kdtree_distance(const double *p1,
const SizeType idx_p2,
std::size_t size) const
{
ValueType s = 0.0;
for (std::size_t i = 0; i < size; i++) {
const ValueType d = p1[i] - m_points[idx_p2][i];
s += d*d;
}
return s;
}
// Returns the dim'th component of the idx'th point in the class:
inline ValueType kdtree_get_pt(const SizeType idx, std::size_t dim) const
{
return m_points[idx][dim];
}
// Optional bounding-box computation: return false to default to a standard bbox computation loop.
// Return true if the BBOX was already computed by the class and returned in "bb" so it can be avoided to redo it again.
// Look at bb.size() to find out the expected dimensionality (e.g. 2 or 3 for point clouds)
template <class BBOX>
bool kdtree_get_bbox(BBOX &bb) const
{
return false;
}
/** @} */
/**
* strategy and nrOfPoints are ignored
* == evalAtPoint_simple
**/
template <class DOFVector>
Value_t<DOFVector> evalAtPoint(DOFVector const& vec, PointType const& x)
{
Value_t<DOFVector> value = 0;
auto const& grid = vec.getFeSpace().gridView().grid();
DataType elementSeed;
if (getNearestElement(x, elementSeed)) {
auto element = grid.entity(elementSeed);
auto lambda = element.geometry().local(x);
value = localEval(vec, element, lambda);
}
return value;
}
private:
/// evaluate DOFVector at barycentric coordinates \p lambda in element
template <class FE, class Value, class Element, class LocalCoords>
Value localEval(DOFVector<FE, Value> const& vec,
Element const& element,
LocalCoords const& lambda)
{
localView.bind(element);
auto const& localFiniteElem = localView.tree().finiteElement();
auto const& localBasis = localFiniteElem.localBasis();
const std::size_t nBasisFct = localFiniteElem.size();
std::vector<Dune::FieldVector<double,1> > shapeValues;
localBasis.evaluateFunction(lambda, shapeValues);
Value data = 0;
for (std::size_t j = 0; j < nBasisFct; ++j) {
auto const global_idx = localView.index(j);
data += vector[global_idx] * shapeValues[j];
}
return data;
}
bool getNearestElement(PointType const& x, DataType& elementSeed)
{
const std::size_t nrOfPoints = 1;
std::vector<DegreeOfFreedom> ret_indexes(nrOfPoints);
std::vector<double> out_dists_sqr(nrOfPoints);
nanoflann::KNNResultSet<double, DegreeOfFreedom> resultSet(nrOfPoints);
resultSet.init(&ret_indexes[0], &out_dists_sqr[0]);
index->findNeighbors(resultSet, x.begin(), nanoflann::SearchParams());
elementSeed = m_data[ret_indexes[0]];
return true;
}
private:
typename FeSpace::LocalView localView;
}; // end of KDTreeMeshAdaptor
template <class FeSpace>
auto provideKDTree(FeSpace const& feSpace)
{
using KD_Tree = KDTreeMeshAdaptor<FeSpace>;
return make_unique<KD_Tree>(feSpace);
}
} // end namespace AMDiS
DIM: 2
DIM_OF_WORLD: 2
number of elements: 4
number of vertices: 5
element vertices:
0 1 4
1 2 4
2 3 4
3 0 4
element boundaries:
0 0 1
0 0 2
0 0 3
0 0 4
vertex coordinates:
0.0 0.0
1.0 0.0
1.0 1.0
0.0 1.0
0.5 0.5
element neighbours:
1 3 -1
2 0 -1
3 1 -1
0 2 -1
{
"dimension of world": 2,
"elliptMesh": {
"macro file name": "macro.stand.2d",
"global refinements": 10
},
"ellipt": {
"mesh": "elliptMesh",
"dim": 2,
"components": 1,
"polynomial degree[0]": 1,
"solver": "cg",
"solver" : {
"max iteration": 1000,
"tolerance": 1e-8,
"info": 10,
"left precon": "diag"
},
"output": {
"filename": "ellipt.2d",
"ParaView format": 1,
"ParaView mode": 1
}
}
}
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
#define TEST_EQ(expr1, expr2) test_exit( (expr1) == (expr2), "'" , #expr1 , "' != '" , #expr2 , "'" )
#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif
using namespace AMDiS;
namespace {
using TestParam = ProblemStatTraits<DIM, DOW, 1>;
using TestProblem = ProblemStat<TestParam>;
using MeshView = typename TestProblem::MeshView;
using Geometry = typename MeshView::template Codim<0>::Geometry;
using RefElements = Dune::ReferenceElements<typename Geometry::ctype, Geometry::mydimension>;
template <class FeSpace>
class LocalEvaluation
{
public:
LocalEvaluation(FeSpace const& feSpace)
: localView(feSpace.localView())
{}
template <class Element>
void bind(Element const& element)
{
localView.bind(element);
}
/// evaluate DOFVector at barycentric coordinates \p lambda in element that
/// is bind to instance in \ref bind().
template <class ValueType, class LocalCoord>
ValueType eval(DOFVector<FeSpace, ValueType> const& vec, LocalCoord const& lambda)
{
auto const& localBasis = localView.tree().finiteElement().localBasis();
localBasis.evaluateFunction(lambda, shapeValues);
ValueType data = 0;
for (std::size_t j = 0; j < shapeValues.size(); ++j) {
const auto global_idx = localView.index(j);
data += vec[global_idx] * shapeValues[j];
}
return data;
}
private:
typename FeSpace::LocalView localView;
std::vector<Dune::FieldVector<double,1> > shapeValues;
};
class AMDiSTester
{
public:
// create a problemStat
AMDiSTester(std::string name)
: prob(make_shared<TestProblem>(name))
{
prob->initialize(INIT_ALL);
}
void test1() // CreateDOFVector
{
using FeSpace = typename TestProblem::template FeSpace<0>;
// construction of DOFVectors
DOFVector<FeSpace, double> vec1(prob->getFeSpace(0_c), "vec1");
auto vec2 = makeDOFVector(prob->getFeSpace(0_c), "vec2");
// retrieving DOFVectors from problemStat
auto vec3 = prob->getSolution(0_c);
auto vec5 = prob->getSolution<0>();
auto vec6 = prob->getSolution(index_<0>());
auto vec7 = prob->getSolution()->getDOFVector(0_c);
auto vec8 = prob->getSolution()->getDOFVector<0>();
auto vec9 = prob->getSolution()->getDOFVector(index_<0>());
auto& sys_vec = *prob->getSolution();
auto vec10 = sys_vec[0_c];
auto vec11 = sys_vec[index_<0>()];
// construction of SystemVector
using FeSpaces = typename TestProblem::FeSpaces;
SystemVector<FeSpaces, double> sys_vec1(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec2 = makeSystemVector(*prob->getFeSpaces(), prob->getComponentNames());
auto sys_vec3 = makeSystemVector(*prob->getFeSpaces(), "sys_vec");
// retrieving systemVector from problemStat
auto sys_vec4 = *prob->getSolution();
}
void test2() // FillDOFVector
{
auto& vec = prob->getSolution(0_c);
auto const& feSpace = vec.getFeSpace();
// interpolate function to dofvector
vec.interpol([](auto const& x) { return x*x; });
// test whether interpolation is fine
using FeSpace = std::decay_t< decltype(feSpace) >;
LocalEvaluation<FeSpace> evaluator(feSpace);
for (auto const& element : elements(feSpace.gridView())) {
evaluator.bind(element);
auto const& refElement = RefElements::general(element.type());
std::size_t nVertices = element.subEntities(DIM);
for (std::size_t i = 0; i < nVertices; ++i) {
auto pos = refElement.position(i, DIM);
auto data = evaluator.eval(vec, pos);
auto x = element.geometry().global(pos);
TEST_EQ(data, x*x);
}
}
}
protected:
shared_ptr<TestProblem> prob;
};
} // end namespace
int main(int argc, char** argv)
{
AMDiS::init(argc, argv);
auto tester = make_shared<AMDiSTester>("test");
tester->test1();
tester->test2();
AMDiS::finalize();
return 0;
}
#include <amdis/AMDiS.hpp>
#include <amdis/ProblemStat.hpp>
#include <amdis/common/Literals.hpp>
#ifndef DIM
#define DIM 2
#endif
#ifndef DOW
#define DOW 2
#endif
using namespace AMDiS;
namespace {
using TestParam = ProblemStatTraits<DIM, DOW, 1>;
using TestProblem = ProblemStat<TestParam>;
using Op = TestProblem::OperatorType;
class AMDiSTester
{
public:
// create a problemStat
AMDiSTester(std::string name)
: prob(make_shared<TestProblem>(name))
{
prob->initialize(INIT_ALL);
}
void test1() // AddOperators
{
auto const& u = prob->getSolution<0>();
// define operators
Op op1;
op1.addZOT( 1.0 );
op1.addZOT( constant(1.0) );
// value of DOFVector at QPs
op1.addZOT( valueOf(u) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }) );
op1.addZOT( valueOfFunc(u, [](auto const& U) { return U; }, 1) );
// partial derivative of DOFVector
op1.addZOT( derivativeOf(u, 0) );
op1.addZOT( derivativeOf(u, 1) );
// op1.addFOT( gradientOf(u), GRD_PSI );
// evaluate functor at coordinates
op1.addZOT( eval([](auto const& x) { return x[0]; }) );
op1.addZOT( [](auto const& x) { return 1.0; } );
// add operator to system
prob->addMatrixOperator(op1, 0, 0);
prob->addVectorOperator(op1, 0);
// define operators with factor argument
Op op2;
op2.addZOT( valueOf(u, 1.0) );
op2.addZOT( derivativeOf(u, 0, 1.0) );
op2.addZOT( derivativeOf(u, 1, 1.0) );
// add operator to system, with optional factor as pointer
double* factor = new double(1.0);
prob->addMatrixOperator(op2, 0, 0, factor);
prob->addVectorOperator(op2, 0, factor);
}
void test2() // MatrixVector
{
auto const& u = prob->getSolution<0>();
using GradOp = std::decay_t< decltype(gradientOf(u)) >;
using DerivOp = std::decay_t< decltype(derivativeOf(u,0)) >;
using Jacobian = typename GradOp::value_type;
using Derivative = typename DerivOp::value_type;
Jacobian v, w;
double erg = v[0].dot(w[0]);
}
protected:
shared_ptr<TestProblem> prob;
};
} // end namespace
int main(int argc, char* argv[])
{
AMDiS::init(argc, argv);
auto tester = make_shared<AMDiSTester>("test");
tester->test1();