Bembel
AdjointDoubleLayerOperator.hpp
1 // This file is part of Bembel, the higher order C++ boundary element library.
2 //
3 // Copyright (C) 2024 see <http://www.bembel.eu>
4 //
5 // It was written as part of a cooperation of J. Doelz, H. Harbrecht, S. Kurz,
6 // M. Multerer, S. Schoeps, and F. Wolf at Technische Universitaet Darmstadt,
7 // Universitaet Basel, and Universita della Svizzera italiana, Lugano. This
8 // source code is subject to the GNU General Public License version 3 and
9 // provided WITHOUT ANY WARRANTY, see <http://www.bembel.eu> for further
10 // information.
11 //
12 #ifndef BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_
13 #define BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_
14 
15 namespace Bembel {
16 // forward declaration of class LaplaceAdjointDoubleLayerOperator in order to
17 // define traits
18 class LaplaceAdjointDoubleLayerOperator;
19 
20 template <>
22  typedef Eigen::VectorXd EigenType;
23  typedef Eigen::VectorXd::Scalar Scalar;
24  enum {
25  OperatorOrder = 0,
26  Form = DifferentialForm::Discontinuous,
27  NumberOfFMMComponents = 1
28  };
29 };
30 
35  : public LinearOperatorBase<LaplaceAdjointDoubleLayerOperator> {
36  // implementation of the kernel evaluation, which may be based on the
37  // information available from the superSpace
38  public:
40  template <class T>
41  void evaluateIntegrand_impl(
42  const T &super_space, const SurfacePoint &p1, const SurfacePoint &p2,
43  Eigen::Matrix<typename LinearOperatorTraits<
45  Eigen::Dynamic, Eigen::Dynamic> *intval) const {
46  auto polynomial_degree = super_space.get_polynomial_degree();
47  auto polynomial_degree_plus_one_squared =
48  (polynomial_degree + 1) * (polynomial_degree + 1);
49 
50  // get evaluation points on unit square
51  auto s = p1.segment<2>(0);
52  auto t = p2.segment<2>(0);
53 
54  // get quadrature weights
55  auto ws = p1(2);
56  auto wt = p2(2);
57 
58  // get points on geometry and tangential derivatives
59  auto x_f = p1.segment<3>(3);
60  auto x_f_dx = p1.segment<3>(6);
61  auto x_f_dy = p1.segment<3>(9);
62  auto y_f = p2.segment<3>(3);
63  auto y_f_dx = p2.segment<3>(6);
64  auto y_f_dy = p2.segment<3>(9);
65 
66  // compute surface measures from tangential derivatives
67  auto x_n = x_f_dx.cross(x_f_dy);
68  auto y_kappa = y_f_dx.cross(y_f_dy).norm();
69 
70  // integrand without basis functions
71  auto integrand = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa * ws * wt;
72 
73  // multiply basis functions with integrand and add to intval, this is an
74  // efficient implementation of
75  // (*intval) += super_space.BasisInteraction(s, t) * evaluateKernel(x_f,
76  // y_f)
77  // * x_kappa * y_kappa * ws * wt;
78  super_space.addScaledBasisInteraction(intval, integrand, s, t);
79 
80  return;
81  }
82 
83  Eigen::Matrix<double, 1, 1> evaluateFMMInterpolation_impl(
84  const SurfacePoint &p1, const SurfacePoint &p2) const {
85  // get evaluation points on unit square
86  auto s = p1.segment<2>(0);
87  auto t = p2.segment<2>(0);
88 
89  // get points on geometry and tangential derivatives
90  auto x_f = p1.segment<3>(3);
91  auto x_f_dx = p1.segment<3>(6);
92  auto x_f_dy = p1.segment<3>(9);
93  auto y_f = p2.segment<3>(3);
94  auto y_f_dx = p2.segment<3>(6);
95  auto y_f_dy = p2.segment<3>(9);
96 
97  // compute surface measure in x and unnormalized normal in y from tangential
98  // derivatives
99  auto x_n = x_f_dx.cross(x_f_dy);
100  auto y_kappa = y_f_dx.cross(y_f_dy).norm();
101 
102  // interpolation
103  Eigen::Matrix<double, 1, 1> intval;
104  intval(0) = evaluateKernelGrad(x_f, x_n, y_f) * y_kappa;
105 
106  return intval;
107  }
108 
112  double evaluateKernelGrad(const Eigen::Vector3d &x,
113  const Eigen::Vector3d &x_n,
114  const Eigen::Vector3d &y) const {
115  auto c = x - y;
116  auto r = c.norm();
117  auto r3 = r * r * r;
118  return -c.dot(x_n) / 4. / BEMBEL_PI / r3;
119  }
120 };
121 
122 } // namespace Bembel
123 #endif // BEMBEL_SRC_LAPLACE_ADJOINTDOUBLELAYEROPERATOR_HPP_
double evaluateKernelGrad(const Eigen::Vector3d &x, const Eigen::Vector3d &x_n, const Eigen::Vector3d &y) const
Gradient of fundamental solution of Laplace problem.
Eigen::Matrix< double, 12, 1 > SurfacePoint
typedef of SurfacePoint
Routines for the evalutation of pointwise errors.
Definition: AnsatzSpace.hpp:14
linear operator base class. this serves as a common interface for existing linear operators.
struct containing specifications on the linear operator has to be specialized or derived for any part...