diff --git a/.gitignore b/.gitignore index d786923ed..1de8e59d2 100644 --- a/.gitignore +++ b/.gitignore @@ -7,4 +7,7 @@ eigen3 .vscode # Important for a nice gh-pages branch node_modules/ -package-lock.json \ No newline at end of file +package-lock.json +# dev +vtune +.vscode diff --git a/Bembel/LinearForm b/Bembel/LinearForm index 836b657f0..8e7229216 100644 --- a/Bembel/LinearForm +++ b/Bembel/LinearForm @@ -24,6 +24,7 @@ #include "src/LinearForm/LinearForm.hpp" #include "src/LinearForm/DirichletTrace.hpp" +#include "src/LinearForm/DirichletTraceOnReferenceDomain.hpp" #include "src/LinearForm/NeumannTrace.hpp" #include "src/LinearForm/DiscreteLinearForm.hpp" #include "src/LinearForm/RotatedTangentialTrace.hpp" diff --git a/Bembel/Spline b/Bembel/Spline index 674c68fba..d286149ec 100644 --- a/Bembel/Spline +++ b/Bembel/Spline @@ -29,6 +29,7 @@ #include #include "src/util/Constants.hpp" +#include "src/util/Macros.hpp" #include "src/LinearOperator/DifferentialFormEnum.hpp" #include "src/Spline/DeBoor.hpp" diff --git a/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp b/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp index 9f6f7316d..116e2ae19 100644 --- a/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp +++ b/Bembel/src/AnsatzSpace/FunctionEvaluator.hpp @@ -100,7 +100,7 @@ class FunctionEvaluator { SurfacePoint sp; ansatz_space_.get_superspace().get_geometry()[patch].updateSurfacePoint( - &sp, ref_point, 1, element.mapToReferenceElement(ref_point)); + &sp, patch, ref_point, 1, element.mapToReferenceElement(ref_point)); return evaluate(element, sp); } Eigen::Matrix< diff --git a/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp b/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp index a32bf2122..2d81f85c4 100644 --- a/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp +++ b/Bembel/src/AnsatzSpace/FunctionEvaluatorEval.hpp @@ -29,8 +29,7 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - return coeff.transpose() * super_space.basis(s) / element.get_h(); + return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h(); } }; @@ -47,15 +46,13 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - auto h = element.get_h(); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); + double h = element.get_h(); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - tangential_coefficients = coeff.transpose() * super_space.basis(s); - return (x_f_dx * tangential_coefficients(0) + - x_f_dy * tangential_coefficients(1)) / + tangential_coefficients = + coeff.transpose() * super_space.basis(p.get_xi()); + return (p.get_f_dx() * tangential_coefficients(0) + + p.get_f_dy() * tangential_coefficients(1)) / h; } @@ -67,14 +64,13 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - auto h = element.get_h(); + double h = element.get_h(); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - phiPhiVec_dx = super_space.basisDx(s); + phiPhiVec_dx = super_space.basisDx(p.get_xi()); Eigen::Matrix::Scalar, Eigen::Dynamic, 1> - phiPhiVec_dy = super_space.basisDy(s); + phiPhiVec_dy = super_space.basisDy(p.get_xi()); return (phiPhiVec_dx.dot(coeff.col(0)) + phiPhiVec_dy.dot(coeff.col(1))) / h / h; } @@ -93,8 +89,7 @@ struct FunctionEvaluatorEval { Scalar, Eigen::Dynamic, getFunctionSpaceVectorDimension()> &coeff) const { - auto s = p.segment<2>(0); - return coeff.transpose() * super_space.basis(s) / element.get_h(); + return coeff.transpose() * super_space.basis(p.get_xi()) / element.get_h(); } }; } // namespace Bembel diff --git a/Bembel/src/AnsatzSpace/Projector.hpp b/Bembel/src/AnsatzSpace/Projector.hpp index 666cff142..4d657b204 100644 --- a/Bembel/src/AnsatzSpace/Projector.hpp +++ b/Bembel/src/AnsatzSpace/Projector.hpp @@ -215,18 +215,16 @@ inline _proj_info makeLocalProjectionTriplets( // Here, we suddenly use the degree for basisevaluation, i.e., // maximal_polynomial_degree-1. This is confusing, but correct and tested. { - double vals_y[Constants::MaxP + 1]; - double vals_x[Constants::MaxP + 1]; for (int iy = 0; iy < maximal_polynomial_degree; ++iy) { - Bembel::Basis::ShapeFunctionHandler::evalBasis( - maximal_polynomial_degree - 1, vals_y, mask[iy]); + Eigen::VectorXd vals_y = Bembel::Basis::ShapeFunctionHandler::evalBasis( + maximal_polynomial_degree - 1, mask[iy]); for (int ix = 0; ix < maximal_polynomial_degree; ++ix) { - Bembel::Basis::ShapeFunctionHandler::evalBasis( - maximal_polynomial_degree - 1, vals_x, mask[ix]); + Eigen::VectorXd vals_x = Bembel::Basis::ShapeFunctionHandler::evalBasis( + maximal_polynomial_degree - 1, mask[ix]); for (int jy = 0; jy < maximal_polynomial_degree; ++jy) { for (int jx = 0; jx < maximal_polynomial_degree; ++jx) { system(iy * masksize + ix, jy * maximal_polynomial_degree + jx) = - vals_x[jx] * vals_y[jy]; + vals_x(jx) * vals_y(jy); } } } diff --git a/Bembel/src/AnsatzSpace/SuperSpace.hpp b/Bembel/src/AnsatzSpace/SuperSpace.hpp index 824da2cbb..fe87462f7 100644 --- a/Bembel/src/AnsatzSpace/SuperSpace.hpp +++ b/Bembel/src/AnsatzSpace/SuperSpace.hpp @@ -15,8 +15,8 @@ namespace Bembel { * \ingroup AnsatzSpace * \brief The superspace manages local polynomial bases on each element of the * mesh and provides an interface to evaluate them. - * - * + * + * */ template struct SuperSpace { @@ -26,7 +26,7 @@ struct SuperSpace { ////////////////////////////////////////////////////////////////////////////// /** * \brief Default constructor for the SuperSpace class. - * + * * This constructor creates a SuperSpace object with default parameters. */ SuperSpace() {} @@ -160,7 +160,8 @@ struct SuperSpace { void map2surface(const ElementTreeNode& e, const Eigen::Vector2d& xi, double w, SurfacePoint* surf_pt) const { Eigen::Vector2d st = e.llc_ + e.get_h() * xi; - mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, st, w, xi); + mesh_->get_geometry()[e.patch_].updateSurfacePoint(surf_pt, e.patch_, st, w, + xi); return; } ////////////////////////////////////////////////////////////////////////////// @@ -196,23 +197,20 @@ struct SuperSpace { void addScaledSurfaceCurlInteraction( Eigen::Matrix* intval, Scalar w, const SurfacePoint& p1, const SurfacePoint& p2) const { - // surface measures - double kappa1 = p1.segment<3>(6).cross(p1.segment<3>(9)).norm(); - double kappa2 = p2.segment<3>(6).cross(p2.segment<3>(9)).norm(); // compute basis functions's surface curl. Each column of s_curl is a basis // function's surface curl at point s. Eigen::Matrix s_curl = - (1.0 / kappa1) * - (-p1.segment<3>(6) * basisDy(p1.segment<2>(0)).transpose() + - p1.segment<3>(9) * basisDx(p1.segment<2>(0)).transpose()); + (1.0 / p1.get_surface_measure()) * + (-p1.get_f_dx() * basisDy(p1.get_xi()).transpose() + + p1.get_f_dy() * basisDx(p1.get_xi()).transpose()); Eigen::Matrix t_curl = - (1.0 / kappa2) * - (-p2.segment<3>(6) * basisDy(p2.segment<2>(0)).transpose() + - p2.segment<3>(9) * basisDx(p2.segment<2>(0)).transpose()); + (1.0 / p2.get_surface_measure()) * + (-p2.get_f_dx() * basisDy(p2.get_xi()).transpose() + + p2.get_f_dy() * basisDx(p2.get_xi()).transpose()); // inner product of surface curls of any two basis functions for (int j = 0; j < polynomial_degree_plus_one_squared; ++j) for (int i = 0; i < polynomial_degree_plus_one_squared; ++i) - (*intval)(j * polynomial_degree_plus_one_squared + i) += + (*intval)(j* polynomial_degree_plus_one_squared + i) += w * s_curl.col(i).dot(t_curl.col(j)); } diff --git a/Bembel/src/Geometry/Patch.hpp b/Bembel/src/Geometry/Patch.hpp index 4c0a121a8..e163a7aec 100644 --- a/Bembel/src/Geometry/Patch.hpp +++ b/Bembel/src/Geometry/Patch.hpp @@ -140,7 +140,6 @@ class Patch { const double tpbasisval = xbasis[i] * ybasis[j]; const int accs = 4 * (numy * (polynomial_degree_x_ * x_location + i) + polynomial_degree_y_ * y_location + j); -#pragma omp simd for (int k = 0; k < 4; k++) tmp[k] += data_[accs + k] * tpbasisval; } } @@ -202,7 +201,6 @@ class Patch { // Here I add up the values of the basis functions in the dc // basis -#pragma omp simd for (int k = 0; k < 4; k++) { tmp[k] += data_[accs + k] * tpbasisval; tmpDx[k] += data_[accs + k] * tpbasisvalDx; @@ -216,13 +214,8 @@ class Patch { delete[] xbasisD; delete[] ybasisD; - Eigen::Matrix out; - - // Eigen::Vector3d out; - double bot = 1. / (tmp[3] * tmp[3]); - -#pragma omp simd + Eigen::Matrix out; for (int k = 0; k < 3; k++) { out(k, 0) = (tmpDx[k] * tmp[3] - tmp[k] * tmpDx[3]) * bot; out(k, 1) = (tmpDy[k] * tmp[3] - tmp[k] * tmpDy[3]) * bot; @@ -278,21 +271,10 @@ class Patch { return evalNormal(Eigen::Vector2d(x, y)); } - // - /** - * \brief Updates the surface point and returns the physical point and the - * derivatives there. - * - * This is a combination of eval und evalJacobian, to avoid duplication of - * work. - * - * \param srf_pt Pointer to the SurfacePoint which gets updated. - * \param ref_pt Point in reference domain with respect to the patch. - * \param w quadrature weight. - * \param xi Point in reference domain with respect to the element. - */ - void updateSurfacePoint(SurfacePoint *srf_pt, - const Eigen::Vector2d &ref_pt, double w, + // This is a combination of eval und evalJacobian, to avoid duplication of + // work. See SurfacePoint.hpp + void updateSurfacePoint(SurfacePoint *srf_pt, const int patch, + const Eigen::Vector2d &ref_pt, const double w, const Eigen::Vector2d &xi) const { const int x_location = Spl::FindLocationInKnotVector(ref_pt(0), unique_knots_x_); @@ -304,11 +286,11 @@ class Patch { const double scaledy = Spl::Rescale(ref_pt(1), unique_knots_y_[y_location], unique_knots_y_[y_location + 1]); - // TODO(Felix) Do not use variable-length arrays in accordance to Google - // Style guide - double *buffer = - new double[2 * (polynomial_degree_x_ + polynomial_degree_y_) + 12]; - for (int i = 0; i < 12; ++i) buffer[i] = 0; + // allocating memory is expensive, use the buffer in SurfacePoint + srf_pt->allocate_buffer(2 * (polynomial_degree_x_ + polynomial_degree_y_) + + 12); + double *buffer = srf_pt->get_buffer(); + std::memset(buffer, 0, 12 * sizeof(double)); double *tmp = buffer; double *tmpDx = tmp + 4; @@ -348,19 +330,19 @@ class Patch { const double bot = 1. / tmp[3]; const double botsqr = bot * bot; - (*srf_pt)(0) = xi(0); - (*srf_pt)(1) = xi(1); - (*srf_pt)(2) = w; - (*srf_pt)(3) = tmp[0] * bot; - (*srf_pt)(4) = tmp[1] * bot; - (*srf_pt)(5) = tmp[2] * bot; - (*srf_pt)(6) = (tmpDx[0] * tmp[3] - tmp[0] * tmpDx[3]) * botsqr; - (*srf_pt)(7) = (tmpDx[1] * tmp[3] - tmp[1] * tmpDx[3]) * botsqr; - (*srf_pt)(8) = (tmpDx[2] * tmp[3] - tmp[2] * tmpDx[3]) * botsqr; - (*srf_pt)(9) = (tmpDy[0] * tmp[3] - tmp[0] * tmpDy[3]) * botsqr; - (*srf_pt)(10) = (tmpDy[1] * tmp[3] - tmp[1] * tmpDy[3]) * botsqr; - (*srf_pt)(11) = (tmpDy[2] * tmp[3] - tmp[2] * tmpDy[3]) * botsqr; - delete[] buffer; + Eigen::Vector3d tmpvec(tmp[0], tmp[1], tmp[2]); + Eigen::Vector3d tmpDxvec(tmpDx[0], tmpDx[1], tmpDx[2]); + Eigen::Vector3d tmpDyvec(tmpDy[0], tmpDy[1], tmpDy[2]); + + srf_pt->set_patch(patch); + srf_pt->set_xi(xi); + srf_pt->set_xi_patch(ref_pt); + srf_pt->set_w(w); + srf_pt->set_f(bot * tmpvec); + srf_pt->set_jacobian(botsqr * (Eigen::Matrix() + << tmpDxvec * tmp[3] - tmpvec * tmpDx[3], + tmpDyvec * tmp[3] - tmpvec * tmpDy[3]) + .finished()); return; } diff --git a/Bembel/src/Geometry/SurfacePoint.hpp b/Bembel/src/Geometry/SurfacePoint.hpp index bbd7e0feb..487727b3b 100644 --- a/Bembel/src/Geometry/SurfacePoint.hpp +++ b/Bembel/src/Geometry/SurfacePoint.hpp @@ -11,35 +11,109 @@ // #ifndef BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ #define BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ + #include + /** * \ingroup Geometry - * \brief typedef of SurfacePoint - * - * This typedef is essential for any evaluation of a bilinear form. It provides - * all required geometry information, stored as follows: - * (0) x-coordinate of the evaluation point in the parameter domain [0,1]^2 - * of the current element, i.e we map [0,1]^2->element->surface - * (1) y-coordinate of the evaluation point in the parameter domain [0,1]^2 - * of the current element, i.e we map [0,1]^2->element->surface - * (2) a quadrature weight. Can be left empty if not used as part of a - * quadrature. - * (3) x-coordinate of patch eval in space - * (4) y-coordinate of patch eval in space - * (5) z-coordinate of patch eval in space - * (6) x-component of derivative in x-dir - * (7) y-component of derivative in x-dir - * (8) z-component of derivative in x-dir - * (9) x-component of derivative in y-dir - * (10) y-component of derivative in y-dir - * (11) z-component of derivative in y-dir - * For application of the pull-back to the reference domain, one requires the - * jacobian of any point on the surface. Calling eval and evalJacobian of the - * Patch class introduces work that needs to be done twice. The - * updateSurdacePoint method is specialized and should be used, since it avoids - * redundant work. + * \brief This class stores quadrature and geometry information on quadrature + *points and provides various methods to compute differential geometric + *entities. **/ -typedef Eigen::Matrix SurfacePoint; +class SurfacePoint { + public: + ~SurfacePoint() { + if (buffer_is_allocated()) delete[] buffer_; + } + + // patch + int get_patch() { return patch_; } + const int get_patch() const { return patch_; } + + // quadrature point + Eigen::Vector2d get_xi() { return xi_; } + const Eigen::Vector2d get_xi() const { return xi_; } + + // quadrature point on patch + Eigen::Vector2d get_xi_patch() { return xi_patch_; } + const Eigen::Vector2d get_xi_patch() const { return xi_patch_; } + + // quadrature weight + double get_w() { return w_; } + const double get_w() const { return w_; } + + // geometry parametrization + Eigen::Vector3d get_f() { return f_; } + const Eigen::Vector3d get_f() const { return f_; } + + // geometry parametrization first derivatives + Eigen::Vector3d get_f_dx() { return jacobian_.col(0); } + const Eigen::Vector3d get_f_dx() const { return jacobian_.col(0); } + Eigen::Vector3d get_f_dy() { return jacobian_.col(1); } + const Eigen::Vector3d get_f_dy() const { return jacobian_.col(1); } + Eigen::Matrix get_jacobian() { return jacobian_; } + const Eigen::Matrix get_jacobian() const { return jacobian_; } + + // normal vectors + Eigen::Vector3d get_normal() { + return jacobian_.col(0).cross(jacobian_.col(1)); + } + const Eigen::Vector3d get_normal() const { + return jacobian_.col(0).cross(jacobian_.col(1)); + } + Eigen::Vector3d get_unit_normal() { return get_normal().normalized(); } + const Eigen::Vector3d get_unit_normal() const { + return get_normal().normalized(); + } + + // surface measure + double get_surface_measure() { return get_normal().norm(); } + const double get_surface_measure() const { return get_normal().norm(); } + + // buffer + bool buffer_is_allocated() { return buffer_is_alloced_; } + const bool buffer_is_allocated() const { return buffer_is_alloced_; } + double *get_buffer() { return buffer_; } + const double *get_buffer() const { return buffer_; } + + // setter + void set_patch(const int patch) { patch_ = patch; } + void set_xi(const Eigen::Vector2d &xi) { xi_ = xi; } + void set_xi_patch(const Eigen::Vector2d &xi_patch) { xi_patch_ = xi_patch; } + void set_w(const double w) { w_ = w; } + void set_f(const Eigen::Vector3d &f) { f_ = f; } + void set_jacobian(const Eigen::Matrix &jacobian) { + jacobian_ = jacobian; + } + void allocate_buffer(const int n) { + assert(n > 0); + if (n != buffer_size_) { + if (buffer_is_allocated()) { + delete[] buffer_; + buffer_is_alloced_ = false; + } + // TODO(Felix) Do not use variable-length arrays in accordance to Google + // Style guide + buffer_ = new double[n]; + buffer_size_ = n; + buffer_is_alloced_ = true; + } + } + + private: + int patch_; + Eigen::Vector2d xi_; + Eigen::Vector2d xi_patch_; + double w_; + Eigen::Vector3d f_; + Eigen::Matrix jacobian_; + bool buffer_is_alloced_ = false; + int buffer_size_ = 0; + double *buffer_ = nullptr; +}; + +typedef std::vector> + ElementSurfacePoints; /** * \ingroup Geometry @@ -50,4 +124,3 @@ typedef std::vector> ElementSurfacePoints; #endif // BEMBEL_SRC_GEOMETRY_SURFACEPOINT_HPP_ - diff --git a/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp b/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp index c5d9d13ef..ee39bc2e8 100644 --- a/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/AdjointDoubleLayerOperator.hpp @@ -43,66 +43,24 @@ class HelmholtzAdjointDoubleLayerOperator 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; + auto integrand = + evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); - // 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); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> 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, 1, 1> intval; - intval(0) = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp index 11407cbda..485093b3f 100644 --- a/Bembel/src/Helmholtz/DoubleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/DoubleLayerOperator.hpp @@ -43,66 +43,24 @@ class HelmholtzDoubleLayerOperator 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_n = y_f_dx.cross(y_f_dy); - // integrand without basis functions - auto integrand = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa * ws * wt; + auto integrand = + evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure() * p1.get_w() * p2.get_w(); - // 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); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> 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_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy); - // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/DoubleLayerPotential.hpp b/Bembel/src/Helmholtz/DoubleLayerPotential.hpp index 4d6ebccc6..44a14bcab 100644 --- a/Bembel/src/Helmholtz/DoubleLayerPotential.hpp +++ b/Bembel/src/Helmholtz/DoubleLayerPotential.hpp @@ -41,26 +41,13 @@ class HelmholtzDoubleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute unnormalized normal from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions // dot: adjoint in first variable - auto integrand = evaluateKernelGrad(point, x_f, x_n) * cauchy_value * ws; + auto integrand = evaluateKernelGrad(point, p.get_f(), p.get_normal()) * + cauchy_value * p.get_w(); return integrand; } diff --git a/Bembel/src/Helmholtz/HypersingularOperator.hpp b/Bembel/src/Helmholtz/HypersingularOperator.hpp index a63d6392a..e072ce59b 100644 --- a/Bembel/src/Helmholtz/HypersingularOperator.hpp +++ b/Bembel/src/Helmholtz/HypersingularOperator.hpp @@ -41,46 +41,27 @@ class HelmholtzHypersingularOperator void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::MatrixXcd *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 - Eigen::Vector2d s = p1.segment<2>(0); - Eigen::Vector2d t = p2.segment<2>(0); - - // get quadrature weights - double ws = p1(2); - double wt = p2(2); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p1.segment<3>(3); - Eigen::Vector3d x_f_dx = p1.segment<3>(6); - Eigen::Vector3d x_f_dy = p1.segment<3>(9); - Eigen::Vector3d y_f = p2.segment<3>(3); - Eigen::Vector3d y_f_dx = p2.segment<3>(6); - Eigen::Vector3d y_f_dy = p2.segment<3>(9); - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy); - Eigen::Vector3d y_n = y_f_dx.cross(y_f_dy); + Eigen::Vector3d x_n = p1.get_normal(); + Eigen::Vector3d y_n = p2.get_normal(); // compute h double h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // evaluate kernel - std::complex kernel = evaluateKernel(x_f, y_f); + std::complex kernel = evaluateKernel(p1.get_f(), p2.get_f()); // integrand without basis functions std::complex integrandScalar = - -kernel * x_n.dot(y_n) * wavenumber2_ * ws * wt; + -kernel * x_n.dot(y_n) * wavenumber2_ * p1.get_w() * p2.get_w(); std::complex integrandCurl = - kernel * x_n.norm() * y_n.norm() * ws * wt / h / h; + kernel * x_n.norm() * y_n.norm() * p1.get_w() * p2.get_w() / h / h; // multiply basis functions with integrand and add to intval, this is an // efficient implementation of - super_space.addScaledBasisInteraction(intval, integrandScalar, s, t); + super_space.addScaledBasisInteraction(intval, integrandScalar, p1.get_xi(), + p2.get_xi()); super_space.addScaledSurfaceCurlInteraction(intval, integrandCurl, p1, p2); return; @@ -88,33 +69,18 @@ class HelmholtzHypersingularOperator Eigen::Matrix, 3, 3> evaluateFMMInterpolation_impl( const SurfacePoint &p1, const SurfacePoint &p2) const { - // get evaluation points on unit square - Eigen::Vector2d s = p1.segment<2>(0); - Eigen::Vector2d t = p2.segment<2>(0); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p1.segment<3>(3); - Eigen::Vector3d x_f_dx = p1.segment<3>(6); - Eigen::Vector3d x_f_dy = p1.segment<3>(9); - Eigen::Vector3d y_f = p2.segment<3>(3); - Eigen::Vector3d y_f_dx = p2.segment<3>(6); - Eigen::Vector3d y_f_dy = p2.segment<3>(9); - - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy); - Eigen::Vector3d y_n = y_f_dx.cross(y_f_dy); - // evaluate kernel - std::complex kernel = evaluateKernel(x_f, y_f); + std::complex kernel = evaluateKernel(p1.get_f(), p2.get_f()); // interpolation Eigen::Matrix, 3, 3> intval; intval.setZero(); - intval(0, 0) = -kernel * wavenumber2_ * x_n.dot(y_n); - intval(1, 1) = kernel * x_f_dy.dot(y_f_dy); - intval(1, 2) = -kernel * x_f_dy.dot(y_f_dx); - intval(2, 1) = -kernel * x_f_dx.dot(y_f_dy); - intval(2, 2) = kernel * x_f_dx.dot(y_f_dx); + intval(0, 0) = + -kernel * wavenumber2_ * p1.get_normal().dot(p2.get_normal()); + intval(1, 1) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); + intval(1, 2) = -kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(2, 1) = -kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(2, 2) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); return intval; } diff --git a/Bembel/src/Helmholtz/SingleLayerOperator.hpp b/Bembel/src/Helmholtz/SingleLayerOperator.hpp index ca87da5c6..b45d7ece7 100644 --- a/Bembel/src/Helmholtz/SingleLayerOperator.hpp +++ b/Bembel/src/Helmholtz/SingleLayerOperator.hpp @@ -47,64 +47,28 @@ class HelmholtzSingleLayerOperator 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(); - // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + std::complex integrand = + evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_surface_measure() * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); // 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); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); return; } Eigen::Matrix, 1, 1> 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 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(); - // interpolation Eigen::Matrix, 1, 1> intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Helmholtz/SingleLayerPotential.hpp b/Bembel/src/Helmholtz/SingleLayerPotential.hpp index 027e4e5f3..908d7d616 100644 --- a/Bembel/src/Helmholtz/SingleLayerPotential.hpp +++ b/Bembel/src/Helmholtz/SingleLayerPotential.hpp @@ -46,29 +46,13 @@ class HelmholtzSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); - // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; - + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp index 0f528b4a8..d171166c9 100644 --- a/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp +++ b/Bembel/src/HomogenisedLaplace/SingleLayerOperator.hpp @@ -20,7 +20,7 @@ class HomogenisedLaplaceSingleLayerOperator; /** * \brief Specification of the LinerOperatorTraits for the Homogenised Laplace. */ -template<> +template <> struct LinearOperatorTraits { typedef Eigen::VectorXd EigenType; typedef Eigen::VectorXd::Scalar Scalar; @@ -36,8 +36,8 @@ struct LinearOperatorTraits { * \brief This class implements the specification of the integration for the * single layer operator for the homogenised Laplace. */ -class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< - HomogenisedLaplaceSingleLayerOperator> { +class HomogenisedLaplaceSingleLayerOperator + : public LinearOperatorBase { // implementation of the kernel evaluation, which may be based on the // information available from the superSpace public: @@ -47,74 +47,34 @@ class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< */ HomogenisedLaplaceSingleLayerOperator() { this->deg = getDegree(HomogenisedLaplaceSingleLayerOperator::precision); - this->cs = getCoefficients( - HomogenisedLaplaceSingleLayerOperator::precision); + this->cs = + getCoefficients(HomogenisedLaplaceSingleLayerOperator::precision); } - 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(); - + template + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, + Eigen::Matrix::Scalar, + Eigen::Dynamic, Eigen::Dynamic> *intval) const { // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + auto integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w(); - // 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); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); 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 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(); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; - + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } @@ -122,9 +82,9 @@ class HomogenisedLaplaceSingleLayerOperator : public LinearOperatorBase< * \brief Fundamental solution of the Homogenised Laplace problem */ double evaluateKernel(const Eigen::Vector3d &x, - const Eigen::Vector3d &y) const { - return k_mod(x - y) - + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); + const Eigen::Vector3d &y) const { + return k_mod(x - y) + + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); } /** diff --git a/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp b/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp index 983276ff2..b8c70d7bf 100644 --- a/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp +++ b/Bembel/src/HomogenisedLaplace/SingleLayerPotential.hpp @@ -15,13 +15,13 @@ namespace Bembel { // forward declaration of class HomogenisedLaplaceSingleLayerPotential // in order to define traits -template +template class HomogenisedLaplaceSingleLayerPotential; /** * \brief Specification of the PotentialTraits for the Homogenised Laplace. */ -template +template struct PotentialTraits> { typedef Eigen::VectorXd::Scalar Scalar; static constexpr int OutputSpaceDimension = 1; @@ -32,51 +32,41 @@ struct PotentialTraits> { * \brief This class implements the specification of the integration for the * single layer potential for the homogenised Laplace. */ -template -class HomogenisedLaplaceSingleLayerPotential : public PotentialBase< - HomogenisedLaplaceSingleLayerPotential, LinOp> { +template +class HomogenisedLaplaceSingleLayerPotential + : public PotentialBase, + LinOp> { // implementation of the kernel evaluation, which may be based on the // information available from the superSpace public: /** - * \brief Constructs an object initialising the coefficients and the degree - * via the static variable HomogenisedLaplaceSingleLayerOperator::precision. - */ + * \brief Constructs an object initialising the coefficients and the degree + * via the static variable HomogenisedLaplaceSingleLayerOperator::precision. + */ HomogenisedLaplaceSingleLayerPotential() { - this->deg = getDegree( - HomogenisedLaplaceSingleLayerOperator::getPrecision()); - this->cs = getCoefficients( - HomogenisedLaplaceSingleLayerOperator::getPrecision()); + this->deg = + getDegree(HomogenisedLaplaceSingleLayerOperator::getPrecision()); + this->cs = + getCoefficients(HomogenisedLaplaceSingleLayerOperator::getPrecision()); } Eigen::Matrix< typename PotentialReturnScalar< - typename LinearOperatorTraits::Scalar, - double>::Scalar, 1, 1> evaluateIntegrand_impl( - const FunctionEvaluator &fun_ev, const ElementTreeNode &element, - const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment < 2 > (0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment < 3 > (3); - auto x_f_dx = p.segment < 3 > (6); - auto x_f_dy = p.segment < 3 > (9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - + typename LinearOperatorTraits::Scalar, double>::Scalar, + 1, 1> + evaluateIntegrand_impl(const FunctionEvaluator &fun_ev, + const ElementTreeNode &element, + const Eigen::Vector3d &point, + const SurfacePoint &p) const { // evaluate kernel - auto kernel = evaluateKernel(point, x_f); + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } @@ -85,9 +75,9 @@ class HomogenisedLaplaceSingleLayerPotential : public PotentialBase< * \brief Fundamental solution of the homogenised Laplace problem */ double evaluateKernel(const Eigen::Vector3d &x, - const Eigen::Vector3d &y) const { - return k_mod(x - y) - + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); + const Eigen::Vector3d &y) const { + return k_mod(x - y) + + evaluate_solid_sphericals(x - y, this->cs, this->deg, false); } private: diff --git a/Bembel/src/Identity/IdentityOperatorBase.hpp b/Bembel/src/Identity/IdentityOperatorBase.hpp index ab3b1d930..c486cc038 100644 --- a/Bembel/src/Identity/IdentityOperatorBase.hpp +++ b/Bembel/src/Identity/IdentityOperatorBase.hpp @@ -28,25 +28,10 @@ class IdentityOperatorBase : public LocalOperatorBase { void evaluateIntegrand_impl(const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::MatrixXd *intval) const { - // get evaluation points on unit square - const auto s = p1.segment<2>(0); - - // get quadrature weights - const auto ws = p1(2); - - // get points on geometry and tangential derivatives - const auto &x_f = p1.segment<3>(3); - const auto &x_f_dx = p1.segment<3>(6); - const auto &x_f_dy = p1.segment<3>(9); - - // compute surface measures from tangential derivatives - const auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // integrand without basis functions - const auto integrand = x_kappa * ws; - - super_space.addScaledBasisInteraction(intval, integrand, s, s); - + const auto integrand = p1.get_surface_measure() * p1.get_w(); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p1.get_xi()); return; } }; diff --git a/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp index 332496778..0901dd225 100644 --- a/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp +++ b/Bembel/src/Laplace/AdjointDoubleLayerOperator.hpp @@ -43,65 +43,25 @@ class LaplaceAdjointDoubleLayerOperator 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; + auto integrand = + evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure() * p1.get_w() * p2.get_w(); + + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); - // 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; + intval(0) = evaluateKernelGrad(p1.get_f(), p1.get_normal(), p2.get_f()) * + p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/DoubleLayerOperator.hpp b/Bembel/src/Laplace/DoubleLayerOperator.hpp index 95260fb5f..97cbef8c4 100644 --- a/Bembel/src/Laplace/DoubleLayerOperator.hpp +++ b/Bembel/src/Laplace/DoubleLayerOperator.hpp @@ -43,68 +43,24 @@ class LaplaceDoubleLayerOperator 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(); - auto y_n = y_f_dx.cross(y_f_dy).normalized(); - // integrand without basis functions auto integrand = - evaluateKernelGrad(x_f, y_f, y_n) * y_kappa * x_kappa * ws * wt; + evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure() * p1.get_w() * p2.get_w(); - // 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); + // multiply basis functions with integrand and add to intval + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); 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_kappa = x_f_dx.cross(x_f_dy).norm(); - auto y_n = y_f_dx.cross(y_f_dy); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernelGrad(x_f, y_f, y_n) * x_kappa; - + intval(0) = evaluateKernelGrad(p1.get_f(), p2.get_f(), p2.get_normal()) * + p1.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/DoubleLayerPotential.hpp b/Bembel/src/Laplace/DoubleLayerPotential.hpp index c4d08cda3..a5d210679 100644 --- a/Bembel/src/Laplace/DoubleLayerPotential.hpp +++ b/Bembel/src/Laplace/DoubleLayerPotential.hpp @@ -41,26 +41,12 @@ class LaplaceDoubleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute unnormalized normal from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy); - // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = - evaluateKernelGrad(point, x_f).dot(x_n) * cauchy_value * ws; + auto integrand = evaluateKernelGrad(point, p.get_f()).dot(p.get_normal()) * + cauchy_value * p.get_w(); return integrand; } diff --git a/Bembel/src/Laplace/HypersingularOperator.hpp b/Bembel/src/Laplace/HypersingularOperator.hpp index 0a5918164..b8ca27eb1 100644 --- a/Bembel/src/Laplace/HypersingularOperator.hpp +++ b/Bembel/src/Laplace/HypersingularOperator.hpp @@ -42,36 +42,13 @@ class LaplaceHypersingularOperator 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; + auto integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w() / h / h; // multiply basis functions with integrand and add to intval, this is an // efficient implementation of @@ -82,28 +59,16 @@ class LaplaceHypersingularOperator 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); + auto kernel = evaluateKernel(p1.get_f(), p2.get_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); + intval(0, 0) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); + intval(0, 1) = -kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(1, 0) = -kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(1, 1) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); return intval; } diff --git a/Bembel/src/Laplace/SingleLayerOperator.hpp b/Bembel/src/Laplace/SingleLayerOperator.hpp index 6a7de23e6..cc47a184f 100644 --- a/Bembel/src/Laplace/SingleLayerOperator.hpp +++ b/Bembel/src/Laplace/SingleLayerOperator.hpp @@ -47,64 +47,28 @@ class LaplaceSingleLayerOperator 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(); - // integrand without basis functions - auto integrand = evaluateKernel(x_f, y_f) * x_kappa * y_kappa * ws * wt; + double integrand = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure() * + p1.get_w() * p2.get_w(); // 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); + super_space.addScaledBasisInteraction(intval, integrand, p1.get_xi(), + p2.get_xi()); 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 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(); - // interpolation Eigen::Matrix intval; - intval(0) = evaluateKernel(x_f, y_f) * x_kappa * y_kappa; + intval(0) = evaluateKernel(p1.get_f(), p2.get_f()) * + p1.get_surface_measure() * p2.get_surface_measure(); return intval; } diff --git a/Bembel/src/Laplace/SingleLayerPotential.hpp b/Bembel/src/Laplace/SingleLayerPotential.hpp index b9978ccb5..7e7b3e950 100644 --- a/Bembel/src/Laplace/SingleLayerPotential.hpp +++ b/Bembel/src/Laplace/SingleLayerPotential.hpp @@ -45,29 +45,13 @@ class LaplaceSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); - // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; - + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp b/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp index c04bca505..dac186379 100644 --- a/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp +++ b/Bembel/src/Laplace/SingleLayerPotentialGradient.hpp @@ -41,28 +41,15 @@ class LaplaceSingleLayerPotentialGradient const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // evaluate kernel - auto kernel = evaluateKernelGrad(point, x_f); + auto kernel = evaluateKernelGrad(point, p.get_f()); // assemble Galerkin solution auto cauchy_value = fun_ev.evaluate(element, p); // integrand without basis functions - auto integrand = kernel * cauchy_value * x_kappa * ws; + auto integrand = + kernel * cauchy_value * p.get_surface_measure() * p.get_w(); return integrand; } diff --git a/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp b/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp index ca96a1a51..613a7111b 100644 --- a/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp +++ b/Bembel/src/LaplaceBeltrami/LaplaceBeltramiOperatorBase.hpp @@ -30,15 +30,9 @@ class LaplaceBeltramiOperatorBase : public LocalOperatorBase { const unsigned int elements_per_direction = (1 << super_space.get_refinement_level()); const double h = 1. / ((double)(elements_per_direction)); - - // compute surface measures from tangential derivatives - double x_kappa = p1.segment<3>(6).cross(p1.segment<3>(9)).norm(); - // integrand without basis functions - double integrand = x_kappa * p1(2) / (h * h); - + double integrand = p1.get_surface_measure() * p1.get_w() / (h * h); super_space.addScaledSurfaceGradientInteraction(intval, integrand, p1, p2); - return; } }; diff --git a/Bembel/src/LinearForm/DirichletTrace.hpp b/Bembel/src/LinearForm/DirichletTrace.hpp index fd6e10e4a..43188c01b 100644 --- a/Bembel/src/LinearForm/DirichletTrace.hpp +++ b/Bembel/src/LinearForm/DirichletTrace.hpp @@ -38,30 +38,11 @@ class DirichletTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *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 = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_kappa = x_f_dx.cross(x_f_dy).norm(); - // integrand without basis functions - auto integrand = function_(x_f) * x_kappa * ws; - + Scalar integrand = + function_(p.get_f()) * p.get_surface_measure() * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); - + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; } diff --git a/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp new file mode 100644 index 000000000..96db5e22a --- /dev/null +++ b/Bembel/src/LinearForm/DirichletTraceOnReferenceDomain.hpp @@ -0,0 +1,56 @@ +// 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_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ +#define BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ + +namespace Bembel { + +template +class DirichletTraceOnReferenceDomain; + +template +struct LinearFormTraits> { + typedef ScalarT Scalar; +}; + +/** + * \ingroup LinearForm + * \brief This class provides an implementation of the Dirichlet trace operator + * and a corresponding method to evaluate the linear form corresponding to the + * right hand side of the system via quadrature. + */ +template +class DirichletTraceOnReferenceDomain + : public LinearFormBase, Scalar> { + public: + DirichletTraceOnReferenceDomain() {} + void set_function( + const std::function &function) { + function_ = function; + } + template + void evaluateIntegrand_impl( + const T &super_space, const SurfacePoint &p, + Eigen::Matrix *intval) const { + // integrand without basis functions + Scalar integrand = function_(p.get_patch(), p.get_xi_patch()) * + p.get_surface_measure() * p.get_w(); + // multiply basis functions with integrand + super_space.addScaledBasis(intval, integrand, p.get_xi()); + return; + } + + private: + std::function function_; +}; +} // namespace Bembel + +#endif // BEMBEL_SRC_LINEARFORM_DIRICHLETTRACEONREFERENCEDOMAIN_HPP_ diff --git a/Bembel/src/LinearForm/NeumannTrace.hpp b/Bembel/src/LinearForm/NeumannTrace.hpp index 1a14a2430..b4936eaf6 100644 --- a/Bembel/src/LinearForm/NeumannTrace.hpp +++ b/Bembel/src/LinearForm/NeumannTrace.hpp @@ -40,30 +40,12 @@ class NeumannTrace : public LinearFormBase, Scalar> { void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *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 = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute (unnormalized) surface normal from tangential derivatives - auto x_f_n = x_f_dx.cross(x_f_dy); - // integrand without basis functions // dot: adjoint in first variable - auto integrand = x_f_n.dot(function_(x_f)) * ws; + auto integrand = p.get_normal().dot(function_(p.get_f())) * p.get_w(); // multiply basis functions with integrand - super_space.addScaledBasis(intval, integrand, s); + super_space.addScaledBasis(intval, integrand, p.get_xi()); return; } diff --git a/Bembel/src/LinearForm/RotatedTangentialTrace.hpp b/Bembel/src/LinearForm/RotatedTangentialTrace.hpp index 69861501b..9ae7e082f 100644 --- a/Bembel/src/LinearForm/RotatedTangentialTrace.hpp +++ b/Bembel/src/LinearForm/RotatedTangentialTrace.hpp @@ -40,45 +40,23 @@ class RotatedTangentialTrace void evaluateIntegrand_impl( const T &super_space, const SurfacePoint &p, Eigen::Matrix *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 = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - auto x_f_dx = p.segment<3>(6); - auto x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - auto x_n = x_f_dx.cross(x_f_dy).normalized(); - // tangential component + quadrature weights // use n x f x n = f-n to avoid troubles with -flto flag in combination // of .cross() - auto fun_x_f = function_(x_f); - auto tangential_component = (fun_x_f - fun_x_f.dot(x_n) * x_n) * ws; + Eigen::Matrix fun_x_f = function_(p.get_f()); + Eigen::Vector3d x_n = p.get_unit_normal(); + Eigen::Matrix tangential_component = + (fun_x_f - fun_x_f.dot(x_n) * x_n) * p.get_w(); // extract tangential component - auto component_x = x_f_dx.dot(tangential_component); - auto component_y = x_f_dy.dot(tangential_component); + Eigen::Matrix components = + p.get_jacobian().transpose() * tangential_component; // evaluate shape functions - auto phiPhiVec = super_space.basis(s); - - // multiply basis functions with integrand - Eigen::Matrix phiPhiMat( - polynomial_degree_plus_one_squared, 2); - phiPhiMat.col(0) = component_x * phiPhiVec; - phiPhiMat.col(1) = component_y * phiPhiVec; + auto phiPhiVec = super_space.basis(p.get_xi()); // compute integrals - (*intval) += phiPhiMat; + (*intval) += phiPhiVec * components.transpose(); return; } diff --git a/Bembel/src/LinearForm/TangentialTrace.hpp b/Bembel/src/LinearForm/TangentialTrace.hpp index fc79f9e6f..4a212ca9e 100644 --- a/Bembel/src/LinearForm/TangentialTrace.hpp +++ b/Bembel/src/LinearForm/TangentialTrace.hpp @@ -43,31 +43,18 @@ class TangentialTrace : public LinearFormBase, Scalar> { int polynomial_degree_plus_one_squared = (polynomial_degree + 1) * (polynomial_degree + 1); - // get evaluation points on unit square - Eigen::Vector2d s = p.segment<2>(0); - - // get quadrature weights - double ws = p(2); - - // get points on geometry and tangential derivatives - Eigen::Vector3d x_f = p.segment<3>(3); - Eigen::Vector3d x_f_dx = p.segment<3>(6); - Eigen::Vector3d x_f_dy = p.segment<3>(9); - - // compute surface measures from tangential derivatives - Eigen::Vector3d x_n = x_f_dx.cross(x_f_dy).normalized(); - // tangential component + quadrature weights - Eigen::Matrix fun_x_f = function_(x_f); - Eigen::Matrix tangential_component = fun_x_f.cross(x_n) * ws; + Eigen::Matrix fun_x_f = function_(p.get_f()); + Eigen::Matrix tangential_component = + fun_x_f.cross(p.get_unit_normal()) * p.get_w(); // extract tangential component - Scalar component_x = x_f_dx.dot(tangential_component); - Scalar component_y = x_f_dy.dot(tangential_component); + Scalar component_x = p.get_f_dx().dot(tangential_component); + Scalar component_y = p.get_f_dy().dot(tangential_component); // evaluate shape functions Eigen::Matrix phiPhiVec = - super_space.basis(s); + super_space.basis(p.get_xi()); // multiply basis functions with integrand Eigen::Matrix phiPhiMat( diff --git a/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp b/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp index f10f39893..5f232c469 100644 --- a/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp +++ b/Bembel/src/LinearOperator/Dummy/DummyOperator.hpp @@ -70,8 +70,8 @@ class DummyOperator : public LinearOperatorBase { const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2, Eigen::Matrix::Scalar, Eigen::Dynamic, Eigen::Dynamic> *intval) const { - (*intval)(0, 0) += - test_func_(p1.segment(3, 2), p2.segment(3, 2)) * p1(2) * p2(2); + (*intval)(0, 0) += test_func_(p1.get_f().head<2>(), p2.get_f().head<2>()) * + p1.get_w() * p2.get_w(); return; } diff --git a/Bembel/src/Maxwell/SingleLayerOperator.hpp b/Bembel/src/Maxwell/SingleLayerOperator.hpp index 15495e744..7de8ed816 100644 --- a/Bembel/src/Maxwell/SingleLayerOperator.hpp +++ b/Bembel/src/Maxwell/SingleLayerOperator.hpp @@ -47,70 +47,40 @@ class MaxwellSingleLayerOperator 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 h = 1. / (1 << super_space.get_refinement_level()); // h = 1 ./ (2^M) // integrand without basis functions - auto kernel_evaluation = evaluateKernel(x_f, y_f) * ws * wt; + auto kernel_evaluation = + evaluateKernel(p1.get_f(), p2.get_f()) * p1.get_w() * p2.get_w(); auto integrand_vector = kernel_evaluation; auto integrand_divergence = -kernel_evaluation / wavenumber2_ / h / h; // vector part: mulitply shape functions with integrand and add to buffer - super_space.addScaledVectorBasisInteraction(intval, integrand_vector, s, t, - x_f_dx, x_f_dy, y_f_dx, y_f_dy); + super_space.addScaledVectorBasisInteraction( + intval, integrand_vector, p1.get_xi(), p2.get_xi(), p1.get_f_dx(), + p1.get_f_dy(), p2.get_f_dx(), p2.get_f_dy()); // divergence part: multiply shape functions with integrand and add to // buffer super_space.addScaledVectorBasisDivergenceInteraction( - intval, integrand_divergence, s, t); + intval, integrand_divergence, p1.get_xi(), p2.get_xi()); return; } Eigen::Matrix, 4, 4> 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); + auto kernel = evaluateKernel(p1.get_f(), p2.get_f()); // interpolation Eigen::Matrix, 4, 4> intval; intval.setZero(); - intval(0, 0) = kernel * x_f_dx.dot(y_f_dx); - intval(0, 2) = kernel * x_f_dx.dot(y_f_dy); - intval(2, 0) = kernel * x_f_dy.dot(y_f_dx); - intval(2, 2) = kernel * x_f_dy.dot(y_f_dy); + intval(0, 0) = kernel * p1.get_f_dx().dot(p2.get_f_dx()); + intval(0, 2) = kernel * p1.get_f_dx().dot(p2.get_f_dy()); + intval(2, 0) = kernel * p1.get_f_dy().dot(p2.get_f_dx()); + intval(2, 2) = kernel * p1.get_f_dy().dot(p2.get_f_dy()); intval(1, 1) = -kernel / wavenumber2_; intval(1, 3) = -kernel / wavenumber2_; intval(3, 1) = -kernel / wavenumber2_; diff --git a/Bembel/src/Maxwell/SingleLayerPotential.hpp b/Bembel/src/Maxwell/SingleLayerPotential.hpp index fbe860425..42bf5071c 100644 --- a/Bembel/src/Maxwell/SingleLayerPotential.hpp +++ b/Bembel/src/Maxwell/SingleLayerPotential.hpp @@ -45,32 +45,20 @@ class MaxwellSingleLayerPotential const ElementTreeNode &element, const Eigen::Vector3d &point, const SurfacePoint &p) const { - // get evaluation points on unit square - auto s = p.segment<2>(0); - - // get quadrature weights - auto ws = p(2); - - // get points on geometry and tangential derivatives - auto x_f = p.segment<3>(3); - // compute surface measures from tangential derivatives auto h = element.get_h(); - // evaluate kernel - auto kernel = evaluateKernel(point, x_f); - auto kernel_gradient = evaluateKernelGrad(point, x_f); - + auto kernel = evaluateKernel(point, p.get_f()); + auto kernel_gradient = evaluateKernelGrad(point, p.get_f()); // assemble Galerkin solution auto scalar_part = fun_ev.evaluate(element, p); auto divergence_part = fun_ev.evaluateDiv(element, p); - // integrand without basis functions, note that the surface measure // disappears for the divergence // auto integrand = kernel * scalar_part * ws; auto integrand = (kernel * scalar_part + 1. / wavenumber2_ * kernel_gradient * divergence_part) * - ws; + p.get_w(); return integrand; } diff --git a/Bembel/src/Spline/Basis.hpp b/Bembel/src/Spline/Basis.hpp index 5e4b88d03..c86116e51 100644 --- a/Bembel/src/Spline/Basis.hpp +++ b/Bembel/src/Spline/Basis.hpp @@ -157,8 +157,8 @@ void Phi_times_Phi_(Eigen::Matrix *c, for (int iy = 0; iy < I; iy++) for (int ix = 0; ix < I; ix++) b[iy * I + ix] = X[ix] * Y[iy]; - for (int i = 0; i < (I * I); i++) - for (int j = 0; j < (I * I); j++) (*c)(i, j) += a[i] * b[j]; + for (int j = 0; j < (I * I); j++) + for (int i = 0; i < (I * I); i++) (*c)(i, j) += a[i] * b[j]; return; } @@ -187,14 +187,14 @@ void Div_Phi_times_Div_Phi_( phiphi_dx_

(&b_dx, 1., eta); phiphi_dy_

(&b_dy, 1., eta); - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i, j) += a_dx[i] * b_dx[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i, j + I2) += a_dx[i] * b_dy[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i + I2, j) += a_dy[i] * b_dx[j]; - for (int i = 0; i < I2; ++i) - for (int j = 0; j < I2; ++j) (*c)(i + I2, j + I2) += a_dy[i] * b_dy[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i, j) += a_dx[i] * b_dx[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i, j + I2) += a_dx[i] * b_dy[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i + I2, j) += a_dy[i] * b_dx[j]; + for (int j = 0; j < I2; ++j) + for (int i = 0; i < I2; ++i) (*c)(i + I2, j + I2) += a_dy[i] * b_dy[j]; return; } diff --git a/Bembel/src/Spline/Bernstein.hpp b/Bembel/src/Spline/Bernstein.hpp index 800ae4ef3..0bc440ed1 100644 --- a/Bembel/src/Spline/Bernstein.hpp +++ b/Bembel/src/Spline/Bernstein.hpp @@ -49,68 +49,38 @@ inline constexpr double Bernstein(double evaluation_point) noexcept { //////////////////////////////////////////////////////////////////////////////// /// Hidden Classes //////////////////////////////////////////////////////////////////////////////// -template +template class HiddenBernsteinClass { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - return in[N] * Bernstein(evaluation_point) + - HiddenBernsteinClass::EvalCoefs(in, evaluation_point); - } - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - return ( - (in[N + 1] - in[N]) * Bernstein(evaluation_point) + - HiddenBernsteinClass::EvalDerCoefs(in, evaluation_point)); - } - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - in[N] += Bernstein(evaluation_point); - HiddenBernsteinClass::EvalBasisPEQ(in, evaluation_point); - return; - } - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - in[N] += (P + 1) * (Bernstein(evaluation_point) - - Bernstein(evaluation_point)); - HiddenBernsteinClass::EvalDerBasisPEQ(in, evaluation_point); - return; - } - static inline void EvalBasis(T *in, double evaluation_point) noexcept { + // think twice when replacing double* by something else. All other attempts + // were slower + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { in[N] = Bernstein(evaluation_point); - HiddenBernsteinClass::EvalBasis(in, evaluation_point); + HiddenBernsteinClass::EvalBasis(in, evaluation_point); return; } - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { + // think twice when replacing double* by something else. All other attempts + // were slower + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { in[N] = (P + 1) * (Bernstein(evaluation_point) - Bernstein(evaluation_point)); - HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); + HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); return; } }; -template -class HiddenBernsteinClass { +template +class HiddenBernsteinClass<0, P> { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - return in[0] * Bernstein<0, P>(evaluation_point); - } - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - return (in[1] - in[0]) * Bernstein<0, P>(evaluation_point); - } - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - in[0] += Bernstein<0, P>(evaluation_point); - return; - } - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - in[0] += (-P - 1) * Bernstein<0, P>(evaluation_point); - return; - } - - static inline void EvalBasis(T *in, double evaluation_point) noexcept { + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { in[0] = Bernstein<0, P>(evaluation_point); return; } - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { in[0] = (-P - 1) * Bernstein<0, P>(evaluation_point); return; } @@ -118,128 +88,43 @@ class HiddenBernsteinClass { // This specialization is needed to get a specialized recursion anchor for the // case P = 0. -template -class HiddenBernsteinClass { +template +class HiddenBernsteinClass<-1, P> { public: - static inline T EvalCoefs(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - assert( - false && - "Pos.A This should not happen. Something is wrong with the recursion"); - }; - static inline T EvalDerCoefs(T *in, double evaluation_point) noexcept { - // P needs to be passed lower to avoid infinite recursion - (void)in; - (void)evaluation_point; - return 0; - }; - static inline void EvalBasis(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; + static inline void EvalBasis(double *in, + const double evaluation_point) noexcept { + BEMBEL_UNUSED_(in); + BEMBEL_UNUSED_(evaluation_point); assert( false && "Pos.C This should not happen. Something is wrong with the recursion"); }; - static inline void EvalDerBasis(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - // P needs to be passed lower to avoid infinite recursion - return; - }; - static inline void EvalBasisPEQ(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; + static inline void EvalDerBasis(double *in, + const double evaluation_point) noexcept { + BEMBEL_UNUSED_(in); + BEMBEL_UNUSED_(evaluation_point); assert( false && "Pos.C This should not happen. Something is wrong with the recursion"); }; - static inline void EvalDerBasisPEQ(T *in, double evaluation_point) noexcept { - (void)in; - (void)evaluation_point; - // P needs to be passed lower to avoid infinite recursion - return; - }; }; //////////////////////////////////////////////////////////////////////////////// /// Evaluation Routines //////////////////////////////////////////////////////////////////////////////// -template -T EvalBernstein(T *in, double evaluation_point) noexcept { - return HiddenBernsteinClass::EvalCoefs(in, evaluation_point); -} - -template -void EvalBernstein(T *in, const std::vector &evaluation_points, - T *out) noexcept { - const int N = evaluation_points.size(); - for (int i = 0; i < N; i++) - out[i] = HiddenBernsteinClass::EvalCoefs(in, evaluation_points[i]); - return; -} - -template -std::vector EvalBernstein( - T *in, const std::vector &evaluation_points) noexcept { - const int N = evaluation_points.size(); - std::vector out(N); - for (int i = 0; i < N; i++) - out[i] = HiddenBernsteinClass::EvalCoefs(in, evaluation_points[i]); - return out; -} -template -void EvalBernsteinBasisPEQ(T *in, double evaluation_point) noexcept { - HiddenBernsteinClass::EvalBasisPEQ(in, evaluation_point); - return; -} - -template -void EvalBernsteinBasis(T *in, double evaluation_point) noexcept { - HiddenBernsteinClass::EvalBasis(in, evaluation_point); - return; -} -//////////////////////////////////////////////////////////////////////////////// -/// Evaluation of the Derivatives -//////////////////////////////////////////////////////////////////////////////// -template -T EvalBernsteinDer(T *in, double evaluation_point) noexcept { - return P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_point); -} - -template -void EvalBernsteinDer(T *in, const std::vector &evaluation_points, - T *out) noexcept { - const int N = evaluation_points.size(); - for (int i = 0; i < N; i++) - out[i] = P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_points[i]); - return; -} - -template -std::vector EvalBernsteinDer( - T *in, const std::vector &evaluation_points) noexcept { - const int N = evaluation_points.size(); - std::vector out(N); - for (int i = 0; i < N; i++) - out[i] = P * HiddenBernsteinClass::EvalDerCoefs( - in, evaluation_points[i]); - return out; -} - -template -void EvalBernsteinDerBasisPEQ(T *in, double evaluation_point) noexcept { - in[P] += P * Bernstein

(evaluation_point); - HiddenBernsteinClass::EvalDerBasisPEQ(in, evaluation_point); - return; +// think twice when replacing double* by something else. All other attempts +// were slower +template +void EvalBernsteinBasis(double *in, const double evaluation_point) noexcept { + HiddenBernsteinClass::EvalBasis(in, evaluation_point); } -template -void EvalBernsteinDerBasis(T *in, double evaluation_point) noexcept { +// think twice when replacing double* by something else. All other attempts +// were slower +template +void EvalBernsteinDerBasis(double *in, const double evaluation_point) noexcept { in[P] = P * Bernstein

(evaluation_point); - HiddenBernsteinClass::EvalDerBasis(in, evaluation_point); + HiddenBernsteinClass

::EvalDerBasis(in, evaluation_point); return; } diff --git a/Bembel/src/Spline/Bezierextraction.hpp b/Bembel/src/Spline/Bezierextraction.hpp index 9fb5086fe..e039ae5b1 100644 --- a/Bembel/src/Spline/Bezierextraction.hpp +++ b/Bembel/src/Spline/Bezierextraction.hpp @@ -53,24 +53,15 @@ Eigen::SparseMatrix MakeProjection(const std::vector &x_knots, Eigen::Matrix tempsolver = Eigen::Matrix::Zero(size_phi, size_phi); - double *coefficients_x = new double[polynomial_degree_x]; - double *coefficients_y = new double[polynomial_degree_y]; - - for (int i = 0; i < polynomial_degree_x; i++) coefficients_x[i] = 0; - - for (int i = 0; i < polynomial_degree_y; i++) coefficients_y[i] = 0; - - auto BezBasis = [](std::vector uniq, int deg, int pos, double pt, - double *coef) { + auto BezBasis = [](std::vector uniq, int deg, int pos, double pt) { std::div_t loc = std::div(pos, deg); if (pt > uniq[loc.quot + 1] || pt < uniq[loc.quot]) { return 0.; } else { - coef[loc.rem] = 1; - double out = Bembel::Basis::ShapeFunctionHandler::evalCoef( - deg - 1, coef, Rescale(pt, uniq[loc.quot], uniq[loc.quot + 1])); - coef[loc.rem] = 0; - return out; + Eigen::VectorXd unit = Eigen::VectorXd::Zero(deg); + unit(loc.rem) = 1.; + return unit.dot(Bembel::Basis::ShapeFunctionHandler::evalBasis( + deg - 1, Rescale(pt, uniq[loc.quot], uniq[loc.quot + 1]))); } }; @@ -81,18 +72,13 @@ Eigen::SparseMatrix MakeProjection(const std::vector &x_knots, for (int y = 0; y < size_phi_y; y++) { for (int x = 0; x < size_phi_x; x++) { tempsolver(y * size_phi_x + x, ix * size_phi_y + iy) = - BezBasis(x_unique_knots, polynomial_degree_x, ix, xpoint[x], - coefficients_x) * - BezBasis(y_unique_knots, polynomial_degree_y, iy, ypoint[y], - coefficients_y); + BezBasis(x_unique_knots, polynomial_degree_x, ix, xpoint[x]) * + BezBasis(y_unique_knots, polynomial_degree_y, iy, ypoint[y]); } } } } - delete[] coefficients_x; - delete[] coefficients_y; - Eigen::FullPivLU> lu_decomp(tempsolver); assert(lu_decomp.rank() == size_phi); /// The inverse matrix, for the solition of the linear system to come. diff --git a/Bembel/src/Spline/Localize.hpp b/Bembel/src/Spline/Localize.hpp index efa0cf49a..275095613 100644 --- a/Bembel/src/Spline/Localize.hpp +++ b/Bembel/src/Spline/Localize.hpp @@ -64,18 +64,14 @@ inline std::vector MakeInterpolationPoints( * \brief returns the coefficients to represent a function in the Bernstein * basis on [0,1]. **/ -inline Eigen::Matrix GetInterpolationMatrix( - int polynomial_degree, const std::vector &mask) { - Eigen::Matrix interpolationMatrix(polynomial_degree + 1, - polynomial_degree + 1); - - double val[Constants::MaxP + 1]; - for (int j = 0; j < polynomial_degree + 1; j++) { - Bembel::Basis::ShapeFunctionHandler::evalBasis(polynomial_degree, val, - mask[j]); - for (int i = 0; i < polynomial_degree + 1; i++) - interpolationMatrix(j, i) = val[i]; - } +inline Eigen::MatrixXd GetInterpolationMatrix(int polynomial_degree, + const std::vector &mask) { + Eigen::MatrixXd interpolationMatrix(polynomial_degree + 1, + polynomial_degree + 1); + + for (int j = 0; j < polynomial_degree + 1; j++) + interpolationMatrix.row(j) = Bembel::Basis::ShapeFunctionHandler::evalBasis( + polynomial_degree, mask[j]); return interpolationMatrix.inverse(); } diff --git a/Bembel/src/Spline/ShapeFunctions.hpp b/Bembel/src/Spline/ShapeFunctions.hpp index abee62cfc..989dff090 100644 --- a/Bembel/src/Spline/ShapeFunctions.hpp +++ b/Bembel/src/Spline/ShapeFunctions.hpp @@ -15,9 +15,6 @@ namespace Bembel { namespace Basis { -using funptr_doubleOut_doubleptrDoubleIn = double (*)(double*, double); -using funptr_voidOut_doubleptrDoubleIn = void (*)(double*, double); - /** * \ingroup Spline * \brief These routines implement a template recursion that allows to choose a @@ -27,79 +24,92 @@ using funptr_voidOut_doubleptrDoubleIn = void (*)(double*, double); template class PSpecificShapeFunctionHandler { public: - inline static double evalCoef(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernstein(ar, x) - : PSpecificShapeFunctionHandler

::evalCoef(p, ar, x); - } - inline static double evalDerCoef(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernsteinDer(ar, x) - : PSpecificShapeFunctionHandler

::evalDerCoef(p, ar, x); - } - inline static void evalBasis(int p, double* ar, double x) { - return p == P ? Bembel::Basis::EvalBernsteinBasis(ar, x) + /** + * \brief Evaluates the Bernstein basis of degree p at point x into ar. + */ + inline static void evalBasis(const int p, double* ar, const double x) { + return p == P ? Bembel::Basis::EvalBernsteinBasis

(ar, x) : PSpecificShapeFunctionHandler

::evalBasis(p, ar, x); } - inline static void evalDerBasis(int p, double* ar, double x) { + /** + * \brief Evaluates the Bernstein basis of degree p at point x into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::Matrix evalBasis( + const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalBasis(p, eval.data(), x); + return eval; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p at point + * x into ar. + */ + inline static void evalDerBasis(const int p, double* ar, const double x) { return p == P - ? Bembel::Basis::EvalBernsteinDerBasis(ar, x) + ? Bembel::Basis::EvalBernsteinDerBasis

(ar, x) : PSpecificShapeFunctionHandler

::evalDerBasis(p, ar, x); } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalCoef( - int p) { - return p == P ? &Bembel::Basis::EvalBernstein - : PSpecificShapeFunctionHandler

::ptrEvalCoef(p); - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalDerCoef( - int p) { - return p == P ? &Bembel::Basis::EvalBernsteinDer - : PSpecificShapeFunctionHandler

::ptrEvalDerCoef(p); - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalBasis(int p) { - return p == P ? &Bembel::Basis::EvalBernsteinBasis - : PSpecificShapeFunctionHandler

::ptrEvalBasis(p); - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalDerBasis( - int p) { - return p == P ? &Bembel::Basis::EvalBernsteinDerBasis - : PSpecificShapeFunctionHandler

::ptrEvalDerBasis(p); - } - inline static constexpr bool checkP(int p) { - static_assert(P > 0, "Polynomial degree must be larger than zero"); - return p <= Constants::MaxP; + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::Matrix evalDerBasis( + const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalDerBasis(p, eval.data(), x); + return eval; } }; template <> class PSpecificShapeFunctionHandler<0> { public: - inline static double evalCoef(int p, double* ar, double x) { - return Bembel::Basis::EvalBernstein(ar, x); - } - inline static double evalDerCoef(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinDer(ar, x); - } - inline static void evalBasis(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinBasis(ar, x); - } - inline static void evalDerBasis(int p, double* ar, double x) { - return Bembel::Basis::EvalBernsteinDerBasis(ar, x); - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalCoef( - int p) { - return &Bembel::Basis::EvalBernstein; - } - inline static constexpr funptr_doubleOut_doubleptrDoubleIn ptrEvalDerCoef( - int p) { - return &Bembel::Basis::EvalBernsteinDer; - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalBasis(int p) { - return &Bembel::Basis::EvalBernsteinBasis; - } - inline static constexpr funptr_voidOut_doubleptrDoubleIn ptrEvalDerBasis( - int p) { - return &Bembel::Basis::EvalBernsteinDerBasis; + /** + * \brief Evaluates the Bernstein basis of degree p at point x into ar. + */ + inline static void evalBasis(const int p, double* ar, const double x) { + Bembel::Basis::EvalBernsteinBasis<0>(ar, x); + return; + } + /** + * \brief Evaluates the Bernstein basis of degree p at point x into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::VectorXd evalBasis(const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalBasis(p, eval.data(), x); + return eval; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p at point + * x into ar. + */ + inline static void evalDerBasis(const int p, double* ar, const double x) { + Bembel::Basis::EvalBernsteinDerBasis<0>(ar, x); + return; + } + /** + * \brief Evaluates the derivative of the Bernstein basis of degree p into an + * Eigen::VectorXd. + * + * \attention For performance critical applications use the pointer version of + * this function to avoid slow memory allocation. + */ + inline static Eigen::VectorXd evalDerBasis(const int p, const double x) { + Eigen::VectorXd eval(p + 1); + evalDerBasis(p, eval.data(), x); + return eval; } - inline static constexpr bool checkP(int p) { return Constants::MaxP >= 0; } }; using ShapeFunctionHandler = PSpecificShapeFunctionHandler; diff --git a/Bembel/src/util/surfaceL2error.hpp b/Bembel/src/util/surfaceL2error.hpp index 5a140a36f..e0f24c21e 100644 --- a/Bembel/src/util/surfaceL2error.hpp +++ b/Bembel/src/util/surfaceL2error.hpp @@ -25,21 +25,18 @@ double surfaceL2error(const AnsatzSpace &ansatz_space, GaussSquare GS; auto Q = GS[deg]; SurfacePoint qp; + const auto longvec = (ansatz_space.get_transformation_matrix() * vec).eval(); const auto &super_space = ansatz_space.get_superspace(); const ElementTree &et = super_space.get_mesh().get_element_tree(); for (auto element = et.cpbegin(); element != et.cpend(); ++element) { for (auto i = 0; i < Q.w_.size(); ++i) { super_space.map2surface(*element, Q.xi_.col(i), Q.w_(i), &qp); - // get points on geometry and tangential derivatives - const auto &x_f = qp.segment<3>(3); - const auto &x_f_dx = qp.segment<3>(6); - const auto &x_f_dy = qp.segment<3>(9); - const auto &normal = x_f_dx.cross(x_f_dy); - // compute surface measures from tangential derivatives - Scalar x_kappa = normal.norm(); + // integrand without basis functions const Scalar val = fun_val.evaluate(*element, qp)(0); - retval += x_kappa * Q.w_(i) * element->get_h() * element->get_h() * - (functor(x_f) - val) * (functor(x_f) - val); + retval += qp.get_surface_measure() * Q.w_(i) * element->get_h() * + element->get_h() * + (functor(qp.get_f()) - val) * + (functor(qp.get_f()) - val); } } return sqrt(retval); diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 20f1b4a07..ae06f892b 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -11,6 +11,7 @@ configure_file(${CMAKE_SOURCE_DIR}/geo/cube_small.dat ${CMAKE_CURRENT_BINARY_DIR set(CIFILES AnsatzSpace BlockClusterTree + DirichletTraceOnReferenceDomain Geometry MassMatrix LaplaceBeltrami diff --git a/examples/DirichletTraceOnReferenceDomain.cpp b/examples/DirichletTraceOnReferenceDomain.cpp new file mode 100644 index 000000000..395de90bb --- /dev/null +++ b/examples/DirichletTraceOnReferenceDomain.cpp @@ -0,0 +1,89 @@ +// 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 "Bembel/src/util/surfaceL2error.hpp" +#include "examples/Error.hpp" + +int main() { + using namespace Bembel; + using namespace Eigen; + + int polynomial_degree_max = 3; + int refinement_level_max = 3; + + std::function fun = [](const Vector3d &in) { + return in(0); + }; + Geometry geometry("sphere.dat"); + 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) { + std::cout << "Degree " << polynomial_degree << " Level " + << refinement_level << std::endl; + // Build ansatz space + AnsatzSpace ansatz_space(geometry, refinement_level, + polynomial_degree); + + // Set up and compute discrete operator + DiscreteLocalOperator disc_op(ansatz_space); + disc_op.compute(); + DiscreteLinearForm, MassMatrixScalarDisc> disc_lf( + ansatz_space); + disc_lf.get_linear_form().set_function(fun); + disc_lf.compute(); + + // Set up solver + SparseLU, COLAMDOrdering> solver; + solver.analyzePattern(disc_op.get_discrete_operator()); + solver.factorize(disc_op.get_discrete_operator()); + VectorXd x = solver.solve(disc_lf.get_discrete_linear_form()); + + // Set up function on reference domain + FunctionEvaluator evaluator(ansatz_space); + evaluator.set_function(x); + std::function fun_disc = + [&](int p, const Vector2d &x) { + return evaluator.evaluateOnPatch(p, x)(0); + }; + + // Set up linear form on reference domain + DiscreteLinearForm, + MassMatrixScalarDisc> + disc_lf2(ansatz_space); + disc_lf2.get_linear_form().set_function(fun_disc); + disc_lf2.compute(); + + // check correctness + assert(std::abs((disc_lf.get_discrete_linear_form() - + disc_lf2.get_discrete_linear_form()) + .norm() / + disc_lf.get_discrete_linear_form().norm()) < 1e-10); + } + + std::cout << std::endl; + } + // The VTKwriter sets up initial geomety information. + std::cout << "\n" << std::string(60, '=') << "\n"; + + return 0; +} diff --git a/tests/test_Bernstein.cpp b/tests/test_Bernstein.cpp index 0e876d41d..5b11a86f7 100644 --- a/tests/test_Bernstein.cpp +++ b/tests/test_Bernstein.cpp @@ -15,25 +15,23 @@ int main() { using namespace Bembel; + using namespace Eigen; // We test the Bernstein with given control points against the deBoor code constexpr int P = Bembel::Constants::MaxP; - double* coefs = new double[P + 1]; - - for (int i = 0; i < P + 1; ++i) { - coefs[i] = (i + 1) / static_cast(P + 1); - } - - Eigen::Map coefs_vector(coefs, 1, P + 1); + VectorXd buffer(P + 1); + VectorXd coefs = + VectorXd::LinSpaced(P + 1, 1, P + 1) / static_cast(P + 1); for (int p = 0; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { - double result1 = Basis::ShapeFunctionHandler::evalCoef(p, coefs, x); + buffer.setZero(); + Basis::ShapeFunctionHandler::evalBasis(p, buffer.data(), x); + double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = - Spl::DeBoor(Eigen::MatrixXd(coefs_vector.leftCols(p + 1)), - Spl::MakeBezierKnotVector(p + 1), v)(0); + double result2 = Spl::DeBoor(MatrixXd(coefs.head(p + 1).transpose()), + Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < Test::Constants::coefficient_accuracy); @@ -43,18 +41,18 @@ int main() { // Now, we do the same for the derivatives for (int p = 1; p <= Bembel::Constants::MaxP; ++p) { for (auto x : Test::Constants::eq_points) { - double result1 = Basis::ShapeFunctionHandler::evalDerCoef(p, coefs, x); + buffer.setZero(); + Basis::ShapeFunctionHandler::evalDerBasis(p, buffer.data(), x); + double result1 = buffer.dot(coefs); std::vector v = {x}; - double result2 = - Spl::DeBoorDer(Eigen::MatrixXd(coefs_vector.leftCols(p + 1)), - Spl::MakeBezierKnotVector(p + 1), v)(0); + double result2 = Spl::DeBoorDer(MatrixXd(coefs.head(p + 1).transpose()), + Spl::MakeBezierKnotVector(p + 1), v)(0); BEMBEL_TEST_IF(std::abs(result1 - result2) < Test::Constants::coefficient_accuracy); } } - delete[] coefs; return 0; } diff --git a/tests/test_SurfacePointUpdate.cpp b/tests/test_SurfacePointUpdate.cpp index 166ce7ff7..6e8c154d0 100644 --- a/tests/test_SurfacePointUpdate.cpp +++ b/tests/test_SurfacePointUpdate.cpp @@ -16,31 +16,25 @@ int main() { using namespace Bembel; + // initialize geometry Test::TestGeometryWriter::writeScreen(); - Bembel::Geometry geometry("test_Screen.dat"); assert(geometry.get_geometry().size() == 1); + // initialize tolerance + double tol = Test::Constants::test_tolerance_geometry; + for (auto x : Test::Constants::eq_points) { for (auto y : Test::Constants::eq_points) { - auto pt = Eigen::Vector2d(x, y); - SurfacePoint srf_pt, srf_pt_ref; - - srf_pt_ref.head(2) = pt; - srf_pt_ref(2) = 3.1415; - srf_pt_ref.segment(3, 3) = geometry.get_geometry()[0].eval(pt); - auto dummy = geometry.get_geometry()[0].evalJacobian(pt); - srf_pt_ref.segment(6, 3) = dummy.col(0); - srf_pt_ref.segment(9, 3) = dummy.col(1); - - auto point = geometry.get_geometry()[0].eval(pt); - auto jacobian = geometry.get_geometry()[0].evalJacobian(pt); - - geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, pt, 3.1415, pt); - - if ((srf_pt - srf_pt_ref).norm() > - Test::Constants::test_tolerance_geometry) - return 1; + Eigen::Vector2d pt = Eigen::Vector2d(x, y); + SurfacePoint srf_pt; + geometry.get_geometry()[0].updateSurfacePoint(&srf_pt, 0, pt, 3.1415, pt); + assert((srf_pt.get_xi() - pt).norm() < tol); + assert(std::abs(srf_pt.get_w() - 3.1415) < tol); + assert((srf_pt.get_f().head<2>() - pt).norm() < tol); + assert(std::abs(srf_pt.get_f()(2) < tol)); + assert((srf_pt.get_f_dx() - Eigen::Vector3d(1., 0., 0.)).norm() < tol); + assert((srf_pt.get_f_dy() - Eigen::Vector3d(0, 1., 0.)).norm() < tol); } } return 0;