From: Eldar Khattatov Date: Thu, 28 Dec 2017 17:37:42 +0000 (-0500) Subject: Fix Hdiv seminorm X-Git-Tag: v9.0.0-rc1~613^2 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=refs%2Fpull%2F5668%2Fhead;p=dealii.git Fix Hdiv seminorm --- diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h index ed0806211f..1e68a740e0 100644 --- a/include/deal.II/numerics/vector_tools.templates.h +++ b/include/deal.II/numerics/vector_tools.templates.h @@ -7173,14 +7173,22 @@ namespace VectorTools case Hdiv_seminorm: for (unsigned int q=0; q= dim, + unsigned int idx = 0; + if (weight!=nullptr) + for (; idx 0) + break; + + Assert (n_components >= idx+dim, ExcMessage ("You can only ask for the Hdiv norm for a finite element " "with at least 'dim' components. In that case, this function " - "will take the divergence of the first 'dim' components.")); + "will find the index of the first non-zero weight and take " + "the divergence of the 'dim' components that follow it.")); + Number sum = 0; // take the trace of the derivatives scaled by the weight and square it - for (unsigned int k=0; k::abs_square(sum) * fe_values.JxW(q); } diff = std::sqrt(diff); diff --git a/tests/vector_tools/integrate_difference_05.cc b/tests/vector_tools/integrate_difference_05.cc new file mode 100644 index 0000000000..c74e9cbff5 --- /dev/null +++ b/tests/vector_tools/integrate_difference_05.cc @@ -0,0 +1,155 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2017 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 at +// the top level of the deal.II distribution. +// +// --------------------------------------------------------------------- + + +// Test integrate_difference with Hdiv_seminorm. Specifically this test +// deals with the case when the ComponentSelectFunction is used and +// the components selected for the Hdiv seminorm computations are not +// the first dim components in the Finite Element. + +#include "../tests.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + +using namespace dealii; + +// First dim components: +// f_x = x^2+y(+z), f_y = x^2+y^2, f_z = z+xy +// div(f) = 2x+2y in 2d +// div(f) = 2x+2y+1 in 3d +// Second dim components: +// g(x,y,z) = 2 * f(x,y,z) +template +class Ref : public Function +{ +public: + Ref() : Function(2*dim) {} + + void vector_value (const Point &p, Vector &values) const + { + switch (dim) + { + case 2: + values[0] = p[0]*p[0]+p[1]; + values[1] = p[0]*p[0]+p[1]*p[1]; + values[2] = 2*(p[0]*p[0]+p[1]); + values[3] = 2*(p[0]*p[0]+p[1]*p[1]); + break; + case 3: + values[0] = p[0]*p[0]+p[1]+p[2]; + values[1] = p[0]*p[0]+p[1]*p[1]; + values[2] = p[2]+p[0]*p[1]; + values[3] = 2*(p[0]*p[0]+p[1]+p[2]); + values[4] = 2*(p[0]*p[0]+p[1]*p[1]); + values[5] = 2*(p[2]+p[0]*p[1]); + break; + default: + Assert(false, ExcNotImplemented()); + } + } +}; + + +template +void test(VectorTools::NormType norm, double value) +{ + Triangulation tria; + GridGenerator::hyper_cube(tria); + tria.refine_global(1); + + FESystem fe(FE_Q(3),dim, + FE_Q(3),dim); + DoFHandler dofh(tria); + dofh.distribute_dofs(fe); + + Vector solution (dofh.n_dofs ()); + VectorTools::interpolate(dofh, Ref(), solution); + + Vector cellwise_errors (tria.n_active_cells()); + + const ComponentSelectFunction mask_2 (std::make_pair(dim,2*dim), 2*dim); + VectorTools::integrate_difference (dofh, + solution, + ZeroFunction(2*dim), + cellwise_errors, + QGauss(5), + norm); + + double error = cellwise_errors.l2_norm(); + + const double difference_1 = std::abs(error-value); + deallog << "computed: " << error + << " expected: " << value + << " difference: " << difference_1 + << std::endl; + Assert(difference_1<1e-10, ExcMessage("Error in integrate_difference, first components")); + + VectorTools::integrate_difference (dofh, + solution, + ZeroFunction(2*dim), + cellwise_errors, + QGauss(5), + norm, + &mask_2); + + error = cellwise_errors.l2_norm(); + const double difference_2 = std::abs(error-2.0*value); + deallog << "computed: " << error + << " expected: " << 2.0*value + << " difference: " << difference_2 + << std::endl; + Assert(difference_2<1e-10, ExcMessage("Error in integrate_difference, second components")); + +} + + +template +void test() +{ + deallog << dim << " dimensions, Hdiv_seminorm:" << std::endl; + double true_value = 0; + switch (dim) + { + case 2: + true_value = std::sqrt(14.0/3.0); + break; + case 3: + true_value = std::sqrt(29.0/3.0); + break; + default: + Assert(false, ExcNotImplemented()); + } + + test(VectorTools::Hdiv_seminorm, true_value); + deallog << "OK" << std::endl; +} + + +int main (int argc, char **argv) +{ + initlog(); + test<2>(); + test<3>(); +} diff --git a/tests/vector_tools/integrate_difference_05.output b/tests/vector_tools/integrate_difference_05.output new file mode 100644 index 0000000000..0b13b736ae --- /dev/null +++ b/tests/vector_tools/integrate_difference_05.output @@ -0,0 +1,9 @@ + +DEAL::2 dimensions, Hdiv_seminorm: +DEAL::computed: 2.16025 expected: 2.16025 difference: 4.44089e-16 +DEAL::computed: 4.32049 expected: 4.32049 difference: 8.88178e-16 +DEAL::OK +DEAL::3 dimensions, Hdiv_seminorm: +DEAL::computed: 3.10913 expected: 3.10913 difference: 1.77636e-15 +DEAL::computed: 6.21825 expected: 6.21825 difference: 3.55271e-15 +DEAL::OK