Skip to content
Snippets Groups Projects
Commit ee358257 authored by Simon Praetorius's avatar Simon Praetorius
Browse files

Extract alglib into solver class

parent feda220f
No related tags found
1 merge request!1Extract alglib into solver class
Pipeline #75918 passed
Stage: test
Hide whitespace changes
Inline Side-by-side
.gitignore
+ 11
− 1
View file @ ee358257
build*/
.cache/
.vscode/
libs/**/
compile_commands.json
\ No newline at end of file
compile_commands.json
*.ini
*.dat
*.csv
*.vtu
*.pvd
*.vtp
callgrind.out.*
\ No newline at end of file
dune/meshdist/closestpoint.hh
+ 123
− 116
View file @ ee358257
......@@ -3,143 +3,150 @@
#if DUNE_MESHDIST_HAVE_ALGLIB
#include <alglib/optimization.h>
#else
#include <dune/geometry/utility/algorithms.hh>
#endif
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/utility/algorithms.hh>
namespace Dune::MeshDist {
#if DUNE_MESHDIST_HAVE_ALGLIB
namespace Impl {
template <class K, int dim>
alglib::real_1d_array toAlglib (const Dune::FieldVector<K,dim>& x)
{
alglib::real_1d_array X;
X.setcontent(dim,x.data());
return X;
}
template <class K, int n>
alglib::real_1d_array toAlglib (const Dune::FieldVector<K,n>& x)
{
alglib::real_1d_array X;
X.setcontent(n,x.data());
return X;
}
template <class K, int dim>
Dune::FieldVector<K,dim> fromAlglib (const alglib::real_1d_array& X)
{
Dune::FieldVector<K,dim> x;
for (int i = 0; i < dim; ++i)
x[i] = X[i];
return x;
}
template <class K, int n>
Dune::FieldVector<K,n> fromAlglib (const alglib::real_1d_array& X)
{
Dune::FieldVector<K,n> x;
for (int i = 0; i < n; ++i)
x[i] = X[i];
return x;
}
} // end namespace Impl
template <class F, class DF, class Range>
struct ClosestPointContext
{
F f;
DF df;
Range y;
};
/**
* \brief Using ALGLIB to solve the least-squares problem `|f(x) - y|^2 -> min`,
* subject to the constraint that x is inside the reference element.
*
* \param[in] f Objective function
* \param[in] df Gradient of the objective function f
* \param[in] y Target value
* \param[inout] x0 Initial guess for the solution and the result
* \param[in] opts Additional parameter to control the method, see \ref `GaussNewtonOptions`
* \result Argument `x0` such that `f(x0) = y` in the sense of least squares error.
*/
template <class F, class DF, class Domain,
class Range = std::invoke_result_t<F, Domain>,
class R = typename Dune::FieldTraits<Domain>::real_type>
void closestPoint (const F& f, const DF& df, Range y, Domain& x0,
Dune::Impl::GaussNewtonOptions<R> opts = {})
template <class R, int dim, int dow>
class ClosestPointSolver
{
using namespace alglib;
real_1d_array X0 = Impl::toAlglib(x0);
minlmstate state;
minlmcreatevj(Domain::dimension, Range::dimension, X0, state);
using Context = ClosestPointContext<F,DF,Range>;
Context ctx{f,df,y};
auto objective_function = [](const real_1d_array& X, real_1d_array& Fi, void* ptr)
template <class F, class DF, class Range>
struct ClosestPointContext
{
Context* ctx = (Context*)(ptr);
auto x = Impl::fromAlglib<R,Domain::dimension>(X);
auto f_x = ctx->f(x);
for (int i = 0; i < Range::dimension; ++i)
Fi[i] = Dune::power(f_x[i] - ctx->y[i], 2);
F f;
DF df;
Range y;
};
auto objective_jac = [](const real_1d_array& X, real_1d_array& Fi, real_2d_array& Jij, void* ptr)
{
Context* ctx = (Context*)(ptr);
auto x = Impl::fromAlglib<R,Domain::dimension>(X);
auto f_x = ctx->f(x);
auto df_x = ctx->df(x);
for (int i = 0; i < Range::dimension; ++i)
Fi[i] = Dune::power(f_x[i] - ctx->y[i], 2);
// box constraints: x_i in [0,1]
alglib::real_1d_array bndL_, bndU_;
for (int i = 0; i < Range::dimension; ++i)
for (int j = 0; j < Domain::dimension; ++j)
Jij[i][j] = 2*(f_x[i] - ctx->y[i])*df_x[j][i];
};
// linear constraint
alglib::real_2d_array C_; // sum_i x_i <= 1
alglib::integer_1d_array CT_; // 1: >=, -1: <=
// The following only works for simplices
real_2d_array C; // x >= 0 and sum_i x_i <= 1
C.setlength(Domain::dimension+1, Domain::dimension+1);
for (int i = 0; i < Domain::dimension+1; ++i)
for (int j = 0; j < Domain::dimension+1; ++j)
C[i][j] = R(i == j || i == Domain::dimension);
integer_1d_array CT; // 1: >=, -1: <=
CT.setlength(Domain::dimension+1);
for (int j = 0; j < Domain::dimension+1; ++j)
CT[j] = j < Domain::dimension ? 1 : -1;
minlmsetlc(state, C, CT, Domain::dimension+1);
minlmsetcond(state, opts.absTol, opts.maxIt);
Dune::FieldVector<R,Domain::dimension> s(1);
minlmsetscale(state, Impl::toAlglib(s));
minlmoptimize(state, objective_function, objective_jac, nullptr, &ctx);
minlmreport rep;
minlmresults(state, X0, rep);
x0 = Impl::fromAlglib<R, Domain::dimension>(X0);
}
#else // DUNE_MESHDIST_HAVE_ALGLIB
/**
* \brief Using a simple Gauss-Newton algorithm to solve the least-squares problem
* `|f(x) - y|^2 -> min`, without enforcing the constraint that x is inside the reference element.
*
* \param[in] f Objective function
* \param[in] df Gradient of the objective function f
* \param[in] y Target value
* \param[inout] x0 Initial guess for the solution and the result
* \param[in] opts Additional parameter to control the method, see \ref `GaussNewtonOptions`
* \result Argument `x0` such that `f(x0) = y` in the sense of least squares error.
*/
template <class F, class DF, class Domain,
class Range = std::invoke_result_t<F, Domain>,
class R = typename Dune::FieldTraits<Domain>::real_type>
void closestPoint (const F& f, const DF& df, Range y, Domain& x0,
Dune::Impl::GaussNewtonOptions<R> opts = {})
{
Dune::Impl::gaussNewton(f, df, y, x0, opts);
}
alglib::real_1d_array scaling_;
mutable alglib::minlmstate state_;
#endif // DUNE_MESHDIST_HAVE_ALGLIB
public:
/**
* \brief Construct an initialize a Levenberg-Marquardt solver to find the closest-point
*
* \param[in] gt The geometry type of the domain
* \param[in] pt Additional parameter to control the method, see \ref `GaussNewtonOptions`
*/
explicit ClosestPointSolver (const GeometryType& gt, const Dune::ParameterTree& pt = {})
{
bndL_.setlength(dim);
bndU_.setlength(dim);
for (int j = 0; j < dim; ++j) {
bndL_[j] = 0.0;
bndU_[j] = 1.0;
}
C_.setlength(1, dim+1);
for (int i = 0; i < dim+1; ++i)
C_[0][i] = 1.0;
CT_.setlength(1);
CT_[0] = -1;
scaling_.setlength(dim);
for (int i = 0; i < dim; ++i)
scaling_[i] = 1.0;
alglib::real_1d_array X0;
X0.setlength(dim);
for (int i = 0; i < dim; ++i)
X0[i] = 1.0/3.0;
alglib::minlmcreatevj(dim, dow, X0, state_);
alglib::minlmsetbc(state_, bndL_, bndU_);
if (gt.isSimplex() && dim > 1)
alglib::minlmsetlc(state_, C_, CT_, 1);
alglib::minlmsetcond(state_, pt.get<R>("abs_tol", 1.e-4), pt.get<int>("max_it", 100));
alglib::minlmsetscale(state_, scaling_);
}
/**
* \brief Using ALGLIB to solve the least-squares problem `|f(x) - y|^2 -> min`,
* subject to the constraint that x is inside the reference element.
*
* \param[in] f Objective function
* \param[in] df Gradient of the objective function f
* \param[in] y Target value
* \param[inout] x0 Initial guess for the solution and the result
* \result Argument `x0` such that `f(x0)` is closest to `y` in the sense of least squares error.
*/
template <class F, class DF>
void operator() (const F& f, const DF& df, const FieldVector<R,dow>& y, FieldVector<R,dim>& x0) const
{
alglib::real_1d_array X0 = Impl::toAlglib(x0);
alglib::minlmrestartfrom(state_, X0);
using Context = ClosestPointContext<F,DF,FieldVector<R,dow>>;
Context ctx{f,df,y};
auto objective_function = [](const alglib::real_1d_array& X, alglib::real_1d_array& Fi, void* ptr)
{
Context* ctx = (Context*)(ptr);
auto x = Impl::fromAlglib<R,dim>(X);
auto f_x = ctx->f(x);
for (int i = 0; i < dow; ++i)
Fi[i] = Dune::power(f_x[i] - ctx->y[i], 2);
};
auto objective_jac = [](const alglib::real_1d_array& X, alglib::real_1d_array& Fi, alglib::real_2d_array& Jij, void* ptr)
{
Context* ctx = (Context*)(ptr);
auto x = Impl::fromAlglib<R,dim>(X);
auto f_x = ctx->f(x);
auto df_x = ctx->df(x);
for (int i = 0; i < dow; ++i)
Fi[i] = Dune::power(f_x[i] - ctx->y[i], 2);
for (int i = 0; i < dow; ++i)
for (int j = 0; j < dim; ++j)
Jij[i][j] = 2*(f_x[i] - ctx->y[i])*df_x[j][i];
};
alglib::minlmoptimize(state_, objective_function, objective_jac, nullptr, &ctx);
alglib::minlmreport rep;
alglib::minlmresults(state_, X0, rep);
x0 = Impl::fromAlglib<R,dim>(X0);
}
};
} // end namespace Dune::MeshDist
......
dune/meshdist/entitysearch.hh
+ 10
− 2
View file @ ee358257
......@@ -8,6 +8,7 @@
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/timer.hh>
#include <dune/geometry/mappedgeometry.hh>
#include <dune/geometry/quadraturerules.hh>
......@@ -73,11 +74,15 @@ class EntitySearch
number[i] = std::max<unsigned int>(1,(unsigned int)(diam/box_size));
}
#ifndef NDEBUG
std::cout << "h in [" << h_min << "," << h_max << "]" << std::endl;
#endif
// TODO: Find better way to increase box size
auto extra = (upper-lower)/10;
#ifndef NDEBUG
std::cout << "time(backGrid): " << t.elapsed() << std::endl;
#endif
return BackGrid{lower-extra, upper+extra, number};
}
......@@ -101,12 +106,15 @@ public:
entitySeeds_.at(index).emplace(indexSet.index(e),e.seed());
}
}
#ifndef NDEBUG
std::cout << "time(build entitysearch): " << t.elapsed() << std::endl;
#endif
}
template <class K>
auto elements (const Dune::FieldVector<K,dimension>& x, int boxRadius = 1) const
auto elements (const Dune::FieldVector<K,dimension>& x, const Dune::ParameterTree& pt = {}) const
{
int boxRadius = pt.get<int>("radius", 1);
auto mi = backGrid_.find(x);
auto neighbors = backGrid_.neighbors(mi,boxRadius);
while (boxesEmpty(neighbors))
......
dune/meshdist/entitysearch2.hh
+ 13
− 6
View file @ ee358257
......@@ -8,6 +8,7 @@
#include <vector>
#include <dune/common/fvector.hh>
#include <dune/common/parametertree.hh>
#include <dune/common/timer.hh>
#include <dune/geometry/mappedgeometry.hh>
#include <dune/geometry/quadraturerules.hh>
......@@ -53,8 +54,9 @@ class EntitySearch2
// TODO: Find better way to increase box size
auto extra = (upper-lower)/10;
#ifndef NDEBUG
std::cout << "time(build boundingBox): " << t.elapsed() << std::endl;
#endif
return std::pair{lower-extra, upper+extra};
}
......@@ -72,8 +74,9 @@ class EntitySearch2
auto refElem = referenceElement(e);
centers[indexSet.index(e)] = localX(refElem.position(0,0));
}
#ifndef NDEBUG
std::cout << "time(build centerList): " << t.elapsed() << std::endl;
#endif
return centers;
}
......@@ -87,8 +90,9 @@ class EntitySearch2
{
seeds[indexSet.index(e)] = e.seed();
}
#ifndef NDEBUG
std::cout << "time(build seedlist): " << t.elapsed() << std::endl;
#endif
return seeds;
}
......@@ -103,11 +107,12 @@ public:
{}
template <class K>
auto elements (const Dune::FieldVector<K,dimension>& x, int numElements = 64) const
auto elements (const Dune::FieldVector<K,dimension>& x, const Dune::ParameterTree& pt = {}) const
{
int numElements = pt.get<int>("num_elements", 16);
std::vector<std::size_t> outIndices(numElements);
std::vector<F> outDistancesSq(numElements);
kdtree_.query(x, numElements, outIndices, outDistancesSq);
outDistancesSq_.resize(numElements);
kdtree_.query(x, numElements, outIndices, outDistancesSq_);
return Dune::transformedRangeView(std::move(outIndices), [&](const auto& index)
{
......@@ -120,6 +125,8 @@ private:
std::vector<GlobalCoordinate> points_;
std::vector<EntitySeed> entitySeeds_;
KDTree<F,dimension> kdtree_;
mutable std::vector<F> outDistancesSq_;
};
} // end namespace Dune::MeshDist
......
dune/meshdist/hausdorffdistance.hh
+ 55
− 42
View file @ ee358257
......@@ -7,6 +7,7 @@
#include <limits>
#include <type_traits>
#include <dune/common/parametertree.hh>
#include <dune/common/ftraits.hh>
#include <dune/geometry/mappedgeometry.hh>
#include <dune/geometry/quadraturerules.hh>
......@@ -39,7 +40,8 @@ quadratureRule (const Entity& entity, int quadOrder)
* Return the one-sided Hausdorff distance d(X1,X2) = sup_{x in X1} inf_{y in X2} ||x - y||_2
*/
template <class Distance = LInfDistance<double>, class Parametrization1, class Parametrization2>
std::floating_point auto hausdorffDistance (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto hausdorffDistance (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
auto distance = Distance::init();
using F = typename Distance::value_type;
......@@ -51,46 +53,52 @@ std::floating_point auto hausdorffDistance (const Parametrization1& X1, const Pa
auto entitySearch = Dune::MeshDist::EntitySearch(X2.entitySet().gridView(), X2);
#endif
using EntitySet2 = std::decay_t<decltype(X2.entitySet())>;
using Solver = ClosestPointSolver<F, EntitySet2::LocalCoordinate::dimension, EntitySet2::GlobalCoordinate::dimension>;Solver solver(GeometryTypes::simplex(EntitySet2::LocalCoordinate::dimension), pt.sub("closest_point"));
auto localX1 = localFunction(X1);
auto localX2 = localFunction(X2);
for (auto const& e1 : elements(X1.entitySet().gridView()))
{
localX1.bind(e1);
auto geometry1 = Dune::MappedGeometry(localX1, Dune::Impl::LocalDerivativeGeometry{e1.geometry()});
auto elements2 = entitySearch.elements(geometry1.center(), pt.sub("entity_search"));
for (auto const& qp : Impl::quadratureRule(e1, quadOrder))
{
auto x1 = localX1(qp.position());
const auto& quad = Impl::quadratureRule(e1, pt.get<int>("quad_order", 4));
std::vector<F> minDists(quad.size(), std::numeric_limits<F>::max());
F distX2 = std::numeric_limits<F>::max();
for (auto e2 : entitySearch.elements(x1))
// traverse the closest elements to e1, using the entitysearch
for (auto const& e2 : elements2)
{
localX2.bind(e2);
auto geometry2 = Dune::MappedGeometry(localX2, Dune::Impl::LocalDerivativeGeometry{e2.geometry()});
// Solve a constrained minimization problem to find a point in the referenceElement
auto closestPoint = [&](auto const& x0, auto const& y) {
auto lambda = x0;
solver( // lambda = argmin_l (||y - f(l)||^2)
[&](const auto& l) { return geometry2.global(l); }, // f(x)
[&](const auto& l) { return geometry2.jacobianTransposed(l); }, // grad f(x)
y, lambda);
return lambda;
};
// compute the local coordinate of the closest point in the element for each quad point
for (std::size_t i = 0; i < quad.size(); ++i)
{
localX2.bind(e2);
auto geometry2 = Dune::MappedGeometry(localX2, Dune::Impl::LocalDerivativeGeometry{e2.geometry()});
auto refElem2 = referenceElement(e2);
// TODO: solve a constrained minimization problem with l in referenceElement
auto closestPoint = [&](auto const& y) {
auto lambda = refElem2.position(0,0);
// Dune::Impl::gaussNewton
Dune::MeshDist::closestPoint( // lambda = argmin_l (||y - f(l)||^2)
[&](const auto& l) { return geometry2.global(l); }, // f(x)
[&](const auto& l) { return geometry2.jacobianTransposed(l); }, // grad f(x)
y, lambda, Dune::Impl::GaussNewtonOptions<F>{.maxIt=10, .absTol=1.e-4});
return lambda;
};
// compute the local coordinate of the closest point in the element
auto l2 = closestPoint(x1);
// if (refElem2.checkInside(l2)) {
using std::min;
auto x2 = localX2(l2);
distX2 = min(distX2, Distance::dist(x1, x2));
// }
auto x1 = localX1(quad[i].position());
auto x2 = localX2(closestPoint(quad[i].position(), x1));
using std::min;
minDists[i] = min(minDists[i], Distance::dist(x1, x2));
}
}
auto dS = geometry1.integrationElement(qp.position()) * qp.weight();
Distance::update(distance, distX2, dS);
// finalize the element, either by computing the integral or taking the supremum
for (std::size_t i = 0; i < quad.size(); ++i)
{
auto dS = geometry1.integrationElement(quad[i].position()) * quad[i].weight();
Distance::update(distance, minDists[i], dS);
}
}
......@@ -103,10 +111,11 @@ std::floating_point auto hausdorffDistance (const Parametrization1& X1, const Pa
* with d(A,B) = sup_{x in A} inf_{y in B} ||x - y||_2
*/
template <class Distance = LInfDistance<double>, class Parametrization1, class Parametrization2>
std::floating_point auto symmetricHausdorffDistance (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto symmetricHausdorffDistance (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
using std::max;
return max(hausdorffDistance(X1,X2,quadOrder), hausdorffDistance(X2,X1,quadOrder));
return max(hausdorffDistance(X1,X2,pt), hausdorffDistance(X2,X1,pt));
}
......@@ -133,9 +142,10 @@ std::floating_point auto area (const Parametrization& X)
* Return the mean error between the surfaces X1 and X2: 1/|area(X1)| int_X1 dist(x, X2) dx
*/
template <class F = double, class Parametrization1, class Parametrization2>
std::floating_point auto meanError (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto meanError (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
return hausdorffDistance<L1Distance<F>>(X1,X2,quadOrder) / Impl::area<F>(X1);
return hausdorffDistance<L1Distance<F>>(X1,X2,pt) / Impl::area<F>(X1);
}
......@@ -143,34 +153,37 @@ std::floating_point auto meanError (const Parametrization1& X1, const Parametriz
* Return the mean error between the surfaces X1 and X2: 1/|area(X1)| int_X1 dist(x, X2) dx
*/
template <class F = double, class Parametrization1, class Parametrization2>
std::floating_point auto symmetricMeanError (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto symmetricMeanError (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
using std::max;
return max(
hausdorffDistance<L1Distance<F>>(X1,X2,quadOrder) / Impl::area<F>(X1),
hausdorffDistance<L1Distance<F>>(X2,X1,quadOrder) / Impl::area<F>(X2));
hausdorffDistance<L1Distance<F>>(X1,X2,pt) / Impl::area<F>(X1),
hausdorffDistance<L1Distance<F>>(X2,X1,pt) / Impl::area<F>(X2));
}
/**
* Return the mean error between the surfaces X1 and X2: sqrt( 1/|area(X1)| int_X1 dist(x, X2)^2 dx )
*/
template <class F = double, class Parametrization1, class Parametrization2>
std::floating_point auto meanSquareError (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto meanSquareError (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
using std::sqrt;
hausdorffDistance<L2Distance<F>>(X1,X2,quadOrder) / sqrt(Impl::area<F>(X1));
return hausdorffDistance<L2Distance<F>>(X1,X2,pt) / sqrt(Impl::area<F>(X1));
}
/**
* Return the mean error between the surfaces X1 and X2: sqrt( 1/|area(X1)| int_X1 dist(x, X2)^2 dx )
*/
template <class F = double, class Parametrization1, class Parametrization2>
std::floating_point auto symmetricMeanSquareError (const Parametrization1& X1, const Parametrization2& X2, int quadOrder = 4)
std::floating_point auto symmetricMeanSquareError (const Parametrization1& X1, const Parametrization2& X2,
const Dune::ParameterTree& pt = {})
{
using std::sqrt; using std::max;
return max(
hausdorffDistance<L2Distance<F>>(X1,X2,quadOrder) / sqrt(Impl::area<F>(X1)),
hausdorffDistance<L2Distance<F>>(X2,X1,quadOrder) / sqrt(Impl::area<F>(X2)));
hausdorffDistance<L2Distance<F>>(X1,X2,pt) / sqrt(Impl::area<F>(X1)),
hausdorffDistance<L2Distance<F>>(X2,X1,pt) / sqrt(Impl::area<F>(X2)));
}
......
src/CMakeLists.txt
+ 4
− 0
View file @ ee358257
......@@ -2,6 +2,10 @@ add_executable("ellipse" ellipse.cc)
target_link_dune_default_libraries("ellipse")
target_link_libraries("ellipse" PRIVATE Dune::MeshDist)
add_executable("ellipse_eoc" ellipse_eoc.cc)
target_link_dune_default_libraries("ellipse_eoc")
target_link_libraries("ellipse_eoc" PRIVATE Dune::MeshDist)
if (dune-vtk_FOUND AND dune-foamgrid_FOUND)
add_executable("meshdist" meshdist.cc)
target_link_dune_default_libraries("meshdist")
......
src/ellipse_eoc.cc 0 → 100644
+ 123
− 0
View file @ ee358257
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#include <config.h>
#include <cmath>
#include <iostream>
#include <numbers>
#include <dune/common/parametertree.hh>
#include <dune/common/parametertreeparser.hh>
#include <dune/common/parallel/mpihelper.hh> // An initializer of MPI
#include <dune/common/exceptions.hh> // We use exceptions
#include <dune/foamgrid/foamgrid.hh>
#include <dune/functions/functionspacebases/compositebasis.hh>
#include <dune/functions/functionspacebases/powerbasis.hh>
#include <dune/functions/functionspacebases/lagrangebasis.hh>
#include <dune/functions/functionspacebases/interpolate.hh>
#include <dune/functions/gridfunctions/discreteglobalbasisfunction.hh>
#include <dune/grid/common/gridfactory.hh>
#include <dune/istl/bvector.hh>
#include <dune/meshdist/hausdorffdistance.hh>
using namespace std;
struct Ellipse
{
double a = 1.0, b = 2.0;
Dune::FieldVector<double,2> operator() (const Dune::FieldVector<double,2>& x) const
{
double theta = std::atan2(x[1],x[0]);
return Dune::FieldVector<double,2>({a*std::cos(theta), b*std::sin(theta)});
}
};
int main(int argc, char** argv)
{
Dune::MPIHelper::instance(argc, argv);
Dune::Timer time;
Dune::ParameterTree pt;
std::string inifile = "none.ini";
if (argc > 1) {
inifile = argv[1];
Dune::ParameterTreeParser::readINITree(inifile, pt);
}
using HostGrid = Dune::FoamGrid<1,2>;
using HostGridFactory = Dune::GridFactory<HostGrid>;
HostGridFactory factory;
// ellipse parametrization
double a = 1.0;
double b = 2.0;
unsigned int numPoints = pt.get<int>("grid.num_points", 32);
for (unsigned int i = 0; i < numPoints; ++i)
{
double theta = (2*std::numbers::pi*i)/numPoints;
factory.insertVertex({a*std::cos(theta), b*std::sin(theta)});
}
for (unsigned int i = 0; i < numPoints; ++i)
{
unsigned int j = (i+1)%numPoints;
factory.insertElement(Dune::GeometryTypes::line, {i,j});
}
auto hostGridPtr = factory.createGrid();
int refine = pt.get<int>("grid.refine", 4);
hostGridPtr->globalRefine(refine);
auto ellipse = Ellipse{a,b};
std::vector<double> errsLinf, errsL1, errsL2;
auto vec0 = Dune::BlockVector<double>{};
auto vec1 = Dune::BlockVector<double>{};
for (int level = 1; level <= hostGridPtr->maxLevel(); ++level)
{
std::cout << "level " << level << "/" << hostGridPtr->maxLevel() << std::endl;
using namespace Dune::Functions::BasisFactory;
auto basis0 = makeBasis(hostGridPtr->levelGridView(level-1),
power<HostGrid::dimensionworld>(lagrange<2>(), flatInterleaved()));
auto basis1 = makeBasis(hostGridPtr->levelGridView(level),
power<HostGrid::dimensionworld>(lagrange<2>(), flatInterleaved()));
Dune::Functions::interpolate(basis0,vec0,ellipse);
Dune::Functions::interpolate(basis1,vec1,ellipse);
using Coord = Dune::FieldVector<double,HostGrid::dimensionworld>;
auto X0 = Dune::Functions::makeDiscreteGlobalBasisFunction<Coord>(basis0, vec0);
auto X1 = Dune::Functions::makeDiscreteGlobalBasisFunction<Coord>(basis1, vec1);
time.reset();
errsLinf.push_back(Dune::MeshDist::hausdorffDistance(X1,X0,pt));
std::cout << "time(hausdorffDistance) = " << time.elapsed() << std::endl;
errsL1.push_back(Dune::MeshDist::meanError(X1,X0,pt));
errsL2.push_back(Dune::MeshDist::meanSquareError(X1,X0,pt));
std::cout << "err_Linf = " << errsLinf.back() << std::endl;
}
std::vector<double> eocLinf, eocL1, eocL2;
for (std::size_t i = 1; i < errsLinf.size(); ++i) {
eocLinf.push_back(log(errsLinf[i-1]/errsLinf[i])/log(2));
eocL1.push_back(log(errsL1[i-1]/errsL1[i])/log(2));
eocL2.push_back(log(errsL2[i-1]/errsL2[i])/log(2));
}
{
std::ofstream out("ellipse_eoc.csv", std::ios_base::out);
out << "level;err_Linf;eoc_Linf;err_L1;eoc_L1;err_L2;eoc_L2" << std::endl;
out << 1 << ";" << errsLinf[0] << ";-;" << errsL1[0] << ";-;" << errsL2[0] << ";-" << std::endl;
for (std::size_t i = 1; i < errsLinf.size(); ++i) {
out << i+1 << ";" << errsLinf[i] << ";" << eocLinf[i-1] << ";"
<< errsL1[i] << ";" << eocL1[i-1] << ";"
<< errsL2[i] << ";" << eocL2[i-1] << std::endl;
}
}
}
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