From bda19ffc93733fa346f75cba38c9ef82860e1db0 Mon Sep 17 00:00:00 2001
From: Coulaud <olivier.coulaud@inria.fr>
Date: Wed, 10 Feb 2021 11:35:19 +0100
Subject: [PATCH] Fix bug in the sign of the gradient.
---
experimental/include/scalfmm/matrix_kernels/laplace.hpp | 8 +++++---
1 file changed, 5 insertions(+), 3 deletions(-)
diff --git a/experimental/include/scalfmm/matrix_kernels/laplace.hpp b/experimental/include/scalfmm/matrix_kernels/laplace.hpp
index 38ca60c55..6261df2a3 100644
--- a/experimental/include/scalfmm/matrix_kernels/laplace.hpp
+++ b/experimental/include/scalfmm/matrix_kernels/laplace.hpp
@@ -191,7 +191,7 @@ namespace scalfmm::matrix_kernels::laplace
/// grad_one_over_r\f$ k(x,y) : R -> R^{d}\f$ with \f$d \f$ the space
/// dimension.
///
- /// \f$k(x,y) = \grad | x - y |^{-1}\f$
+ /// \f$k(x,y) = \grad | x - y |^{-1} = -(x-y) | x - y |^{-3}\f$
///
/// scale factor \f$k(ax,ay)= 1/a^2 k(x,y)\f$
template<std::size_t Dim = 3>
@@ -242,7 +242,7 @@ namespace scalfmm::matrix_kernels::laplace
using std::sqrt;
ValueType tmp = ValueType(1.0) / sqrt((((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...));
ValueType r3{pow(tmp, int(3))};
- return std::make_tuple((r3 * (xs.at(Is) - ys.at(Is)))...);
+ return std::make_tuple((r3 * (ys.at(Is) - xs.at(Is)))...);
//-r3 * (xs - ys);
}
@@ -257,6 +257,8 @@ namespace scalfmm::matrix_kernels::laplace
/// val_grad_one_over_r \f$k(x,y) : R^{km} -> R^{kn}\f$ with \f$ km = 1\f$ and \f$kn = d+1\f$ with \f$ d\f$ the
/// space dimension.
///
+ /// \f$k(x,y) = ( | x - y |^{-1}, \grad | x - y |^{-1} ) =(| x - y |^{-1}, -(x-y) | x - y |^{-3}\f)$
+ ///
/// Is a specific kernel used to compute the value of the kernel and its gradient
///
/// scale factor \f$ k(ax,ay)= (1/a, 1/a^2 ... 1/a^2) k(x,y)\f$
@@ -302,7 +304,7 @@ namespace scalfmm::matrix_kernels::laplace
using std::sqrt;
ValueType tmp = ValueType(1.0) / sqrt((((xs.at(Is) - ys.at(Is)) * (xs.at(Is) - ys.at(Is))) + ...));
ValueType r3{pow(tmp, int(3))};
- return std::make_tuple(tmp, (r3 * (xs.at(Is) - ys.at(Is)))...); //-r3 * (xs - ys);
+ return std::make_tuple(tmp, (r3 * (ys.at(Is) - xs.at(Is)))...); //-r3 * (xs - ys);
}
const std::size_t separation_criterion{1};
--
GitLab