From: Luca Heltai Date: Sat, 13 Feb 2021 14:46:22 +0000 (+0100) Subject: Fixed gradient of tensor product cutoff. X-Git-Tag: v9.3.0-rc1~54^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=211a06f03d43678a01064181509c3e3e7ce5ad5f;p=dealii.git Fixed gradient of tensor product cutoff. --- diff --git a/doc/news/changes/minor/20210213LucaHeltai b/doc/news/changes/minor/20210213LucaHeltai new file mode 100644 index 0000000000..bb1e472931 --- /dev/null +++ b/doc/news/changes/minor/20210213LucaHeltai @@ -0,0 +1,3 @@ +Fixed: CutOffFunctionTensorProduct::gradient() was computing the wrong gradient. This is now fixed. +
+(Luca Heltai, 2021/02/13) diff --git a/source/base/function_lib_cutoff.cc b/source/base/function_lib_cutoff.cc index c3fbdb8507..03922d8d46 100644 --- a/source/base/function_lib_cutoff.cc +++ b/source/base/function_lib_cutoff.cc @@ -162,8 +162,13 @@ namespace Functions { Assert(initialized, ExcNotInitialized()); Tensor<1, dim> ret; - for (unsigned int i = 0; i < dim; ++i) - ret[i] = base[i]->gradient(Point<1>(p[i]), component)[0]; + for (unsigned int d = 0; d < dim; ++d) + { + ret[d] = base[d]->gradient(Point<1>(p[d]), component)[0]; + for (unsigned int i = 0; i < dim; ++i) + if (i != d) + ret[d] *= base[i]->value(Point<1>(p[i]), component); + } return ret; } diff --git a/tests/base/function_cutoff_03.cc b/tests/base/function_cutoff_03.cc new file mode 100644 index 0000000000..7c1d3a919e --- /dev/null +++ b/tests/base/function_cutoff_03.cc @@ -0,0 +1,85 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2019 - 2020 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + +// Check that CutOffTensorProductFunction produces the right gradient in dim>1 + +#include +#include +#include + +#include "../tests.h" + +using namespace Functions; + +template +void +test() +{ + Point p; + for (unsigned int i = 0; i < dim; ++i) + p[i] = .5; + + + CutOffFunctionTensorProduct fun( + .5, p, 1, CutOffFunctionBase::no_component, true); + + + fun.template set_base(); + + deallog << "Center: " << fun.get_center() << std::endl + << "Radius: " << fun.get_radius() << std::endl; + + // Check integration of the function, and of its gradient square + QGauss quad(12); + + std::vector> gradients(quad.size()); + std::vector values(quad.size()); + fun.value_list(quad.get_points(), values); + fun.gradient_list(quad.get_points(), gradients); + + // Integral in 1d of grad square: 2.0*pi^2 + // Integral in 2d of grad square: 6.0*pi^2 + // Integral in 3d of grad square: 13.5*pi^2 + + const double spi = numbers::PI * numbers::PI; + + const std::array exact{{2.0 * spi, 6.0 * spi, 13.5 * spi}}; + + double integral = 0.0; + double grad_integral = 0.0; + for (unsigned int i = 0; i < quad.size(); ++i) + { + integral += values[i] * quad.weight(i); + grad_integral += gradients[i] * gradients[i] * quad.weight(i); + } + if (std::abs(integral - 1.0) > 1e-10) + deallog << "Value NOT OK" << std::endl; + else + deallog << "Value OK" << std::endl; + if (std::abs(grad_integral - exact[dim - 1]) > 1e-10) + deallog << "Gradient NOT OK" << std::endl; + else + deallog << "Gradient OK" << std::endl; +} + +int +main() +{ + initlog(1); + + test<1>(); + test<2>(); + test<3>(); +} diff --git a/tests/base/function_cutoff_03.output b/tests/base/function_cutoff_03.output new file mode 100644 index 0000000000..8c8ba15336 --- /dev/null +++ b/tests/base/function_cutoff_03.output @@ -0,0 +1,13 @@ + +DEAL::Center: 0.500000 +DEAL::Radius: 0.500000 +DEAL::Value OK +DEAL::Gradient OK +DEAL::Center: 0.500000 0.500000 +DEAL::Radius: 0.500000 +DEAL::Value OK +DEAL::Gradient OK +DEAL::Center: 0.500000 0.500000 0.500000 +DEAL::Radius: 0.500000 +DEAL::Value OK +DEAL::Gradient OK