From 3dfde86fba768be6a29d2e6e423c45fc534cfc4e Mon Sep 17 00:00:00 2001 From: Maximilian <38978321+mx-nlte@users.noreply.github.com> Date: Tue, 19 Mar 2024 15:17:59 +0100 Subject: [PATCH] Add Missing Integral Operators Laplace with Examples (#94) Hypersingular Operator Double Layer Operator Adjoint Double Layer Operator --- .github/workflows/build-n-deploy.yml | 6 +- Bembel/Laplace | 2 + .../Laplace/AdjointDoubleLayerOperator.hpp | 123 ++++++++++++++ Bembel/src/Laplace/HypersingularOperator.hpp | 153 ++++++++++++++++++ examples/CMakeLists.txt | 3 + examples/Data.hpp | 7 + examples/LaplaceAdjointDoubleLayerH2.cpp | 124 ++++++++++++++ examples/LaplaceDoubleLayerH2.cpp | 116 +++++++++++++ examples/LaplaceHypersingularH2.cpp | 129 +++++++++++++++ 9 files changed, 662 insertions(+), 1 deletion(-) create mode 100644 Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp create mode 100644 Bembel/src/Laplace/HypersingularOperator.hpp create mode 100644 examples/LaplaceAdjointDoubleLayerH2.cpp create mode 100644 examples/LaplaceDoubleLayerH2.cpp create mode 100644 examples/LaplaceHypersingularH2.cpp diff --git a/.github/workflows/build-n-deploy.yml b/.github/workflows/build-n-deploy.yml index f7ec97c2e..1451e743f 100644 --- a/.github/workflows/build-n-deploy.yml +++ b/.github/workflows/build-n-deploy.yml @@ -24,13 +24,17 @@ jobs: - name: Checkout 🛎️ uses: actions/checkout@v4 + - name: Configure system + run: | + sudo apt-get update + - name: Getting Doxygen run: | sudo apt-get install --yes doxygen - name: Install Latex run: | - sudo apt install --yes texlive-latex-extra + sudo apt-get install --yes texlive-base - name: Getting dot run: | diff --git a/Bembel/Laplace b/Bembel/Laplace index 120e1bec7..b532171b0 100644 --- a/Bembel/Laplace +++ b/Bembel/Laplace @@ -27,8 +27,10 @@ #include "Potential" #include "src/Laplace/DoubleLayerOperator.hpp" +#include "src/Laplace/AdjointDoubleLayerOperator.hpp" #include "src/Laplace/DoubleLayerPotential.hpp" #include "src/Laplace/SingleLayerOperator.hpp" +#include "src/Laplace/HypersingularOperator.hpp" #include "src/Laplace/SingleLayerPotential.hpp" #include "src/Laplace/SingleLayerPotentialGradient.hpp" diff --git a/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp new file mode 100644 index 000000000..332496778 --- /dev/null +++ b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp @@ -0,0 +1,123 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. +// +#ifndef BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_ +#define BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_ + +namespace Bembel { +// forward declaration of class LaplaceAdjointDoubleLayerOperator in order to +// define traits +class LaplaceAdjointDoubleLayerOperator; + +template <> +struct LinearOperatorTraits { + typedef Eigen::VectorXd EigenType; + typedef Eigen::VectorXd::Scalar Scalar; + enum { + OperatorOrder = 0, + Form = DifferentialForm::Discontinuous, + NumberOfFMMComponents = 1 + }; +}; + +/** + * \ingroup Laplace + */ +class LaplaceAdjointDoubleLayerOperator + : public LinearOperatorBase { + // implementation of the kernel evaluation, which may be based on the + // information available from the superSpace + public: + LaplaceAdjointDoubleLayerOperator() {} + template + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, + Eigen::Matrix::Scalar, + Eigen::Dynamic, Eigen::Dynamic> *intval) const { + auto polynomial_degree = super_space.get_polynomial_degree(); + auto polynomial_degree_plus_one_squared = + (polynomial_degree + 1) * (polynomial_degree + 1); + + // get evaluation points on unit square + auto s = p1.segment<2>(0); + auto t = p2.segment<2>(0); + + // get quadrature weights + auto ws = p1(2); + auto wt = p2(2); + + // get points on geometry and tangential derivatives + auto x_f = p1.segment<3>(3); + auto x_f_dx = p1.segment<3>(6); + auto x_f_dy = p1.segment<3>(9); + auto y_f = p2.segment<3>(3); + auto y_f_dx = p2.segment<3>(6); + auto y_f_dy = p2.segment<3>(9); + + // compute surface measures from tangential derivatives + auto x_n = x_f_dx.cross(x_f_dy); + auto y_kappa = y_f_dx.cross(y_f_dy).norm(); + + // integrand without basis functions + auto integrand = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa * ws * wt; + + // multiply basis functions with integrand and add to intval, this is an + // efficient implementation of + // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f, + // y_f) + // * x_kappa * y_kappa * ws * wt; + super_space.addScaledBasisInteraction(intval, integrand, s, t); + + return; + } + + Eigen::Matrix evaluateFMMInterpolation_impl( + const SurfacePoint &p1, const SurfacePoint &p2) const { + // get evaluation points on unit square + auto s = p1.segment<2>(0); + auto t = p2.segment<2>(0); + + // get points on geometry and tangential derivatives + auto x_f = p1.segment<3>(3); + auto x_f_dx = p1.segment<3>(6); + auto x_f_dy = p1.segment<3>(9); + auto y_f = p2.segment<3>(3); + auto y_f_dx = p2.segment<3>(6); + auto y_f_dy = p2.segment<3>(9); + + // compute surface measure in x and unnormalized normal in y from tangential + // derivatives + auto x_n = x_f_dx.cross(x_f_dy); + auto y_kappa = y_f_dx.cross(y_f_dy).norm(); + + // interpolation + Eigen::Matrix intval; + intval(0) = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa; + + return intval; + } + + /** + * \brief Gradient of fundamental solution of Laplace problem + */ + double evaluateKernelGrad(const Eigen::Vector3d &x, + const Eigen::Vector3d &x_n, + const Eigen::Vector3d &y) const { + auto c = x - y; + auto r = c.norm(); + auto r3 = r * r * r; + return -c.dot(x_n) / 4. / BEMBEL_PI / r3; + } +}; + +} // namespace Bembel +#endif // BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_ diff --git a/Bembel/src/Laplace/HypersingularOperator.hpp b/Bembel/src/Laplace/HypersingularOperator.hpp new file mode 100644 index 000000000..0a5918164 --- /dev/null +++ b/Bembel/src/Laplace/HypersingularOperator.hpp @@ -0,0 +1,153 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. +#ifndef BEMBEL_SRC_LAPLACE_HYPERSINGULAROPERATOR_HPP_ +#define BEMBEL_SRC_LAPLACE_HYPERSINGULAROPERATOR_HPP_ + +namespace Bembel { +// forward declaration of class LaplaceHypersingularOperator in order to define +// traits +class LaplaceHypersingularOperator; + +template <> +struct LinearOperatorTraits { + typedef Eigen::VectorXd EigenType; + typedef Eigen::VectorXd::Scalar Scalar; + enum { + OperatorOrder = 1, + Form = DifferentialForm::Continuous, + NumberOfFMMComponents = 2 + }; +}; + +/** + * \ingroup Laplace + */ +class LaplaceHypersingularOperator + : public LinearOperatorBase { + // implementation of the kernel evaluation, which may be based on the + // information available from the superSpace + public: + LaplaceHypersingularOperator() {} + template + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, + Eigen::Matrix< + typename LinearOperatorTraits::Scalar, + Eigen::Dynamic, Eigen::Dynamic> *intval) const { + auto polynomial_degree = super_space.get_polynomial_degree(); + auto polynomial_degree_plus_one_squared = + (polynomial_degree + 1) * (polynomial_degree + 1); + + // get evaluation points on unit square + auto s = p1.segment<2>(0); + auto t = p2.segment<2>(0); + + // get quadrature weights + auto ws = p1(2); + auto wt = p2(2); + + // get points on geometry and tangential derivatives + auto x_f = p1.segment<3>(3); + auto x_f_dx = p1.segment<3>(6); + auto x_f_dy = p1.segment<3>(9); + auto y_f = p2.segment<3>(3); + auto y_f_dx = p2.segment<3>(6); + auto y_f_dy = p2.segment<3>(9); + + // compute surface measures from tangential derivatives + auto x_kappa = x_f_dx.cross(x_f_dy).norm(); + auto y_kappa = y_f_dx.cross(y_f_dy).norm(); + + // compute h + auto h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) + + // integrand without basis functions + auto integrand = + evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt / h / h; + + // multiply basis functions with integrand and add to intval, this is an + // efficient implementation of + super_space.addScaledSurfaceCurlInteraction(intval, integrand, p1, p2); + + return; + } + + Eigen::Matrix evaluateFMMInterpolation_impl( + const SurfacePoint &p1, const SurfacePoint &p2) const { + // get evaluation points on unit square + auto s = p1.segment<2>(0); + auto t = p2.segment<2>(0); + + // get points on geometry and tangential derivatives + auto x_f = p1.segment<3>(3); + auto x_f_dx = p1.segment<3>(6); + auto x_f_dy = p1.segment<3>(9); + auto y_f = p2.segment<3>(3); + auto y_f_dx = p2.segment<3>(6); + auto y_f_dy = p2.segment<3>(9); + + // evaluate kernel + auto kernel = evaluateKernel(x_f, y_f); + + // interpolation + Eigen::Matrix intval; + intval.setZero(); + intval(0, 0) = kernel * x_f_dy.dot(y_f_dy); + intval(0, 1) = -kernel * x_f_dy.dot(y_f_dx); + intval(1, 0) = -kernel * x_f_dx.dot(y_f_dy); + intval(1, 1) = kernel * x_f_dx.dot(y_f_dx); + + return intval; + } + + /** + * \brief Fundamental solution of Laplace problem + */ + double evaluateKernel(const Eigen::Vector3d &x, + const Eigen::Vector3d &y) const { + return 1. / 4. / BEMBEL_PI / (x - y).norm(); + } +}; + +/** + * \brief The hypersingular operator requires a special treatment of the + * moment matrices of the FMM due to the involved derivatives on the ansatz + * functions. + */ +template +struct H2Multipole::Moment2D { + static std::vector compute2DMoment( + const SuperSpace &super_space, + const int cluster_level, const int cluster_refinements, + const int number_of_points) { + Eigen::MatrixXd moment_dx = moment2DComputer< + Moment1DDerivative, + Moment1D>( + super_space, cluster_level, cluster_refinements, number_of_points); + Eigen::MatrixXd moment_dy = moment2DComputer< + Moment1D, + Moment1DDerivative>( + super_space, cluster_level, cluster_refinements, number_of_points); + + Eigen::MatrixXd moment(moment_dx.rows() + moment_dy.rows(), + moment_dx.cols()); + moment << moment_dx, moment_dy; + + std::vector vector_of_moments; + vector_of_moments.push_back(moment); + + return vector_of_moments; + } +}; + +} // namespace Bembel +#endif // BEMBEL_SRC_LAPLACE_HYPERSINGULAROPERATOR_HPP_ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index ca5a6d78f..20f1b4a07 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -16,6 +16,9 @@ set(CIFILES LaplaceBeltrami LaplaceSingleLayerFull LaplaceSingleLayerH2 + LaplaceHypersingularH2 + LaplaceDoubleLayerH2 + LaplaceAdjointDoubleLayerH2 LazyEigenSum HelmholtzSingleLayerFull HelmholtzSingleLayerH2 diff --git a/examples/Data.hpp b/examples/Data.hpp index 212d75c3b..1e9e6dd1f 100644 --- a/examples/Data.hpp +++ b/examples/Data.hpp @@ -24,6 +24,13 @@ namespace Data { inline double HarmonicFunction(Eigen::Vector3d in) { return (4 * in(0) * in(0) - 3 * in(1) * in(1) - in(2) * in(2)); } +/* @brief This function implements the gradient of a harmonic function, + * which, in the interior dormain, satisfies the Laplace equation. + * + */ +inline Eigen::Vector3d HarmonicFunctionGrad(Eigen::Vector3d in) { + return Eigen::Vector3d(8 * in(0), -6 * in(1), -2 * in(2)); +} /* @brief This function implements the Helmholtz fundamental solution, * which, if center is placed in the interior domain, satisfies the Helmholtz diff --git a/examples/LaplaceAdjointDoubleLayerH2.cpp b/examples/LaplaceAdjointDoubleLayerH2.cpp new file mode 100644 index 000000000..512c2f49a --- /dev/null +++ b/examples/LaplaceAdjointDoubleLayerH2.cpp @@ -0,0 +1,124 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "examples/Data.hpp" +#include "examples/Error.hpp" +#include "examples/Grids.hpp" + +int main() { + using namespace Bembel; + using namespace Eigen; + + Bembel::IO::Stopwatch sw; + + int polynomial_degree_max = 2; + int refinement_level_max = 3; + + // Load geometry from file "sphere.dat", which must be placed in the same + // directory as the executable + Geometry geometry("sphere.dat"); + + // Define evaluation points for potential field, a tensor product grid of + // 10*10*10 points in [-.25,.25]^3 + MatrixXd gridpoints = Util::makeTensorProductGrid( + VectorXd::LinSpaced(10, -.25, .25), VectorXd::LinSpaced(10, -.25, .25), + VectorXd::LinSpaced(10, -.25, .25)); + + // Define analytical solution using lambda function, in this case a harmonic + // function, see Data.hpp + std::function fun = [](Vector3d in) { + return Data::HarmonicFunction(in); + }; + std::function funGrad = [](Vector3d in) { + return Data::HarmonicFunctionGrad(in); + }; + + std::cout << "\n" << std::string(60, '=') << "\n"; + // Iterate over polynomial degree + for (int polynomial_degree = 0; polynomial_degree < polynomial_degree_max + 1; + ++polynomial_degree) { + VectorXd error(refinement_level_max + 1); + // Iterate over refinement levels + for (int refinement_level = 0; refinement_level < refinement_level_max + 1; + ++refinement_level) { + sw.tic(); + std::cout << "Degree " << polynomial_degree << " Level " + << refinement_level << "\t\t"; + // Build ansatz space + AnsatzSpace ansatz_space_helm( + geometry, refinement_level, polynomial_degree); + AnsatzSpace ansatz_space_mass( + geometry, refinement_level, polynomial_degree); + + // Set up load vector + DiscreteLinearForm, + LaplaceAdjointDoubleLayerOperator> + disc_lf(ansatz_space_helm); + disc_lf.get_linear_form().set_function(funGrad); + disc_lf.compute(); + + // Set up and compute discrete operator + DiscreteOperator, LaplaceAdjointDoubleLayerOperator> + disc_op_double(ansatz_space_helm); + disc_op_double.compute(); + const H2Matrix &AK = disc_op_double.get_discrete_operator(); + + // Assemble mass matrix and system matrix + DiscreteLocalOperator disc_op_mass( + ansatz_space_mass); + disc_op_mass.compute(); + SparseMatrix M = disc_op_mass.get_discrete_operator(); + auto system_matrix = 0.5 * M + AK; // important: do NOT change auto! + + // solve system + GMRES gmres; + gmres.compute(system_matrix); + VectorXd rho = gmres.solve(disc_lf.get_discrete_linear_form()); + + // evaluate potential + DiscretePotential< + LaplaceSingleLayerPotential, + LaplaceAdjointDoubleLayerOperator> + disc_pot(ansatz_space_helm); + disc_pot.set_cauchy_data(rho); + auto pot = disc_pot.evaluate(gridpoints); + + // compute reference, print time, and compute error + VectorXcd pot_ref(gridpoints.rows()); + for (int i = 0; i < gridpoints.rows(); ++i) + pot_ref(i) = fun(gridpoints.row(i)); + error(refinement_level) = (pot - pot_ref).cwiseAbs().maxCoeff(); + std::cout << " time " << std::setprecision(4) << sw.toc() << "s\t\t"; + std::cout << error(refinement_level) << std::endl; + } + + // estimate rate of convergence and check whether it is at least 90% of the + // expected value + assert( + checkRateOfConvergence(error.tail(3), 2 * polynomial_degree + 2, 0.9)); + + std::cout << std::endl; + } + std::cout << std::string(60, '=') << std::endl; + + return 0; +} diff --git a/examples/LaplaceDoubleLayerH2.cpp b/examples/LaplaceDoubleLayerH2.cpp new file mode 100644 index 000000000..08a1889a6 --- /dev/null +++ b/examples/LaplaceDoubleLayerH2.cpp @@ -0,0 +1,116 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "examples/Data.hpp" +#include "examples/Error.hpp" +#include "examples/Grids.hpp" + +int main() { + using namespace Bembel; + using namespace Eigen; + + Bembel::IO::Stopwatch sw; + + int polynomial_degree_max = 2; + int refinement_level_max = 3; + + // Load geometry from file "sphere.dat", which must be placed in the same + // directory as the executable + Geometry geometry("sphere.dat"); + + // Define evaluation points for potential field, a tensor product grid of + // 7*7*7 points in [-.1,.1]^3 + MatrixXd gridpoints = Util::makeTensorProductGrid( + VectorXd::LinSpaced(10, -.25, .25), VectorXd::LinSpaced(10, -.25, .25), + VectorXd::LinSpaced(10, -.25, .25)); + + // Define analytical solution using lambda function, in this case a harmonic + // function, see Data.hpp + std::function fun = [](Vector3d in) { + return Data::HarmonicFunction(in); + }; + + std::cout << "\n" << std::string(60, '=') << "\n"; + // Iterate over polynomial degree + for (int polynomial_degree = 0; polynomial_degree < polynomial_degree_max + 1; + ++polynomial_degree) { + VectorXd error(refinement_level_max + 1); + // Iterate over refinement levels + for (int refinement_level = 0; refinement_level < refinement_level_max + 1; + ++refinement_level) { + sw.tic(); + std::cout << "Degree " << polynomial_degree << " Level " + << refinement_level << "\t\t"; + // Build ansatz space + AnsatzSpace ansatz_space_helm( + geometry, refinement_level, polynomial_degree); + AnsatzSpace ansatz_space_mass( + geometry, refinement_level, polynomial_degree); + + // Set up load vector + DiscreteLinearForm, LaplaceDoubleLayerOperator> + disc_lf(ansatz_space_helm); + disc_lf.get_linear_form().set_function(fun); + disc_lf.compute(); + + // Set up and compute discrete operator + DiscreteOperator, LaplaceDoubleLayerOperator> + disc_op_double(ansatz_space_helm); + disc_op_double.compute(); + const H2Matrix &K = + disc_op_double.get_discrete_operator(); + + // Assemble mass matrix and system matrix + DiscreteLocalOperator disc_op_mass( + ansatz_space_mass); + disc_op_mass.compute(); + SparseMatrix M = disc_op_mass.get_discrete_operator(); + auto system_matrix = -0.5 * M + K; // important: do NOT change auto! + + // solve system + GMRES gmres; + gmres.compute(system_matrix); + VectorXd rho = gmres.solve(disc_lf.get_discrete_linear_form()); + + // evaluate potential + DiscretePotential, + LaplaceDoubleLayerOperator> + disc_pot(ansatz_space_helm); + disc_pot.set_cauchy_data(rho); + auto pot = disc_pot.evaluate(gridpoints); + + error(refinement_level) = maxPointwiseError(pot, gridpoints, fun); + std::cout << " time " << std::setprecision(4) << sw.toc() << "s\t\t"; + std::cout << error(refinement_level) << std::endl; + } + + // estimate rate of convergence and check whether it is at least 90% of the + // expected value + assert( + checkRateOfConvergence(error.tail(3), 2 * polynomial_degree + 2, 0.9)); + + std::cout << std::endl; + } + std::cout << std::string(60, '=') << std::endl; + + return 0; +} diff --git a/examples/LaplaceHypersingularH2.cpp b/examples/LaplaceHypersingularH2.cpp new file mode 100644 index 000000000..1ca68d976 --- /dev/null +++ b/examples/LaplaceHypersingularH2.cpp @@ -0,0 +1,129 @@ +// This file is part of Bembel, the higher order C++ boundary element library. +// +// Copyright (C) 2024 see +// +// It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz, +// M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt, +// Universitaet Basel, and Universita della Svizzera italiana, Lugano. This +// source code is subject to the GNU General Public License version 3 and +// provided WITHOUT ANY WARRANTY, see for further +// information. + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "examples/Data.hpp" +#include "examples/Error.hpp" +#include "examples/Grids.hpp" + +int main() { + using namespace Bembel; + using namespace Eigen; + + Bembel::IO::Stopwatch sw; + + int polynomial_degree_max = 2; + int refinement_level_max = 3; + + // Load geometry from file "sphere.dat", which must be placed in the same + // directory as the executable + Geometry geometry("sphere.dat"); + + // Define evaluation points for potential field, a tensor product grid of + // 7*7*7 points in [-.1,.1]^3 + MatrixXd gridpoints = Util::makeTensorProductGrid( + VectorXd::LinSpaced(10, -.25, .25), VectorXd::LinSpaced(10, -.25, .25), + VectorXd::LinSpaced(10, -.25, .25)); + + // Define analytical solution using lambda function, in this case a harmonic + // function, see Data.hpp + std::function fun = [](Vector3d in) { + return Data::HarmonicFunction(in); + }; + std::function funGrad = [](Vector3d in) { + return Data::HarmonicFunctionGrad(in); + }; + std::function one = [](Vector3d in) { return 1.; }; + + std::cout << "\n" << std::string(60, '=') << "\n"; + // Iterate over polynomial degree + for (int polynomial_degree = 1; polynomial_degree < polynomial_degree_max + 1; + ++polynomial_degree) { + VectorXd error(refinement_level_max + 1); + // Iterate over refinement levels + for (int refinement_level = 0; refinement_level < refinement_level_max + 1; + ++refinement_level) { + sw.tic(); + + std::cout << "Degree " << polynomial_degree << " Level " + << refinement_level << "\t\t"; + // Build ansatz space + AnsatzSpace ansatz_space( + geometry, refinement_level, polynomial_degree); + + // Set up load vector + DiscreteLinearForm, LaplaceHypersingularOperator> + disc_lf(ansatz_space); + disc_lf.get_linear_form().set_function(funGrad); + disc_lf.compute(); + + // Set up one + DiscreteLinearForm, LaplaceHypersingularOperator> + one_lf(ansatz_space); + one_lf.get_linear_form().set_function(one); + one_lf.compute(); + + // Set up and compute discrete operator + DiscreteOperator, LaplaceHypersingularOperator> disc_op( + ansatz_space); + disc_op.compute(); + + // solve system + ConjugateGradient, Lower | Upper, IdentityPreconditioner> + cg; + cg.compute(disc_op.get_discrete_operator()); + VectorXd rho = cg.solve(-disc_lf.get_discrete_linear_form()); + + // evaluate potential + DiscretePotential< + LaplaceDoubleLayerPotential, + LaplaceHypersingularOperator> + disc_pot(ansatz_space); + disc_pot.set_cauchy_data(rho); + VectorXd pot = disc_pot.evaluate(gridpoints); + + // factor out constant + pot -= + 0.5 * VectorXd::Ones(pot.rows()) * (pot.maxCoeff() + pot.minCoeff()); + + // compute reference, print time, and compute error + VectorXcd pot_ref(gridpoints.rows()); + for (int i = 0; i < gridpoints.rows(); ++i) + pot_ref(i) = fun(gridpoints.row(i)); + error(refinement_level) = (pot - pot_ref).cwiseAbs().maxCoeff(); + + // print time and error + std::cout << " time " << std::setprecision(4) << sw.toc() << "s\t\t"; + std::cout << error(refinement_level) << std::endl; + } + + // estimate rate of convergence and check whether it is at least 90% of the + // expected value + assert( + checkRateOfConvergence(error.tail(3), 2 * polynomial_degree + 1, 0.9)); + + std::cout << std::endl; + } + std::cout << "=============================================================" + "==========" + << std::endl; + + return 0; +}