From b217217091c80d0126b7c7f8ffca852fdccb1c24 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Mon, 3 Feb 2014 17:06:24 +0000 Subject: [PATCH] Fix a bug reported by Christoph Heiniger (including a description of how it should be fixed): The computation of curls was wrong. git-svn-id: https://svn.dealii.org/trunk@32394 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/doc/news/changes.h | 6 ++++ deal.II/source/fe/fe_values.cc | 41 +++++++++++++++----------- tests/deal.II/fe_values_view_28.cc | 2 +- tests/deal.II/fe_values_view_28.output | 35 ++++++---------------- 4 files changed, 40 insertions(+), 44 deletions(-) diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h index 56b56bf4dc..be608b1720 100644 --- a/deal.II/doc/news/changes.h +++ b/deal.II/doc/news/changes.h @@ -124,6 +124,12 @@ inconvenience this causes.

Specific improvements

    +
  1. Fixed: Calling FEValuesViews::Vector::get_function_curls() computed + wrong results in some cases (see https://code.google.com/p/dealii/issues/detail?id=182). + This is now fixed. +
    + (Christoph Heiniger, Wolfgang Bangerth, 2014/02/03) +
  2. Added: The class LAPACKFullMatrix now implements interfaces to matrix-matrix multiplication. Also, LAPACKFullMatrix::apply_lu_factorization now also operates on multiple right hand sides in form of another diff --git a/deal.II/source/fe/fe_values.cc b/deal.II/source/fe/fe_values.cc index aaf9fe0a0a..163f03d2b7 100644 --- a/deal.II/source/fe/fe_values.cc +++ b/deal.II/source/fe/fe_values.cc @@ -800,24 +800,22 @@ namespace FEValuesViews const dealii::Tensor<1, spacedim> *shape_gradient_ptr = &shape_gradients[snc][0]; - switch (shape_function_data[shape_function].single_nonzero_component_index) - { - case 0: - { - for (unsigned int q_point = 0; - q_point < n_quadrature_points; ++q_point) - curls[q_point][0] -= value * (*shape_gradient_ptr++)[1]; - - break; - } - - default: - for (unsigned int q_point = 0; - q_point < n_quadrature_points; ++q_point) - curls[q_point][0] += value * (*shape_gradient_ptr)[0]; - } + Assert (shape_function_data[shape_function].single_nonzero_component >= 0, + ExcInternalError()); + // we're in 2d, so the formula for the curl is simple: + if (shape_function_data[shape_function].single_nonzero_component_index == 0) + for (unsigned int q_point = 0; + q_point < n_quadrature_points; ++q_point) + curls[q_point][0] -= value * (*shape_gradient_ptr++)[1]; + else + for (unsigned int q_point = 0; + q_point < n_quadrature_points; ++q_point) + curls[q_point][0] += value * (*shape_gradient_ptr++)[0]; } else + // we have multiple non-zero components in the shape functions. not + // all of them must necessarily be within the 2-component window + // this FEValuesViews::Vector object considers, however. { if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { @@ -887,17 +885,26 @@ namespace FEValuesViews break; } - default: + case 2: + { for (unsigned int q_point = 0; q_point < n_quadrature_points; ++q_point) { curls[q_point][0] += value * (*shape_gradient_ptr)[1]; curls[q_point][1] -= value * (*shape_gradient_ptr++)[0]; } + break; + } + + default: + Assert (false, ExcInternalError()); } } else + // we have multiple non-zero components in the shape functions. not + // all of them must necessarily be within the 3-component window + // this FEValuesViews::Vector object considers, however. { if (shape_function_data[shape_function].is_nonzero_shape_function_component[0]) { diff --git a/tests/deal.II/fe_values_view_28.cc b/tests/deal.II/fe_values_view_28.cc index 74b1736bcd..284011d540 100644 --- a/tests/deal.II/fe_values_view_28.cc +++ b/tests/deal.II/fe_values_view_28.cc @@ -124,7 +124,7 @@ void test_hyper_cube() Triangulation tr; GridGenerator::hyper_cube(tr); - FE_Nedelec fe(3); + FESystem fe(FE_Q(2), dim); test(tr, fe); } diff --git a/tests/deal.II/fe_values_view_28.output b/tests/deal.II/fe_values_view_28.output index 9341da7943..c3c1e10f0b 100644 --- a/tests/deal.II/fe_values_view_28.output +++ b/tests/deal.II/fe_values_view_28.output @@ -1,27 +1,10 @@ -DEAL::FE=FESystem<2>[FE_Q<2>(1)^2] -DEAL:: curls[q]= -2.00 -DEAL:: grads[q]= 2.00 4.00 2.00 4.00 -DEAL:: curls[q]= -2.00 -DEAL:: grads[q]= 2.00 4.00 2.00 4.00 -DEAL:: curls[q]= -2.00 -DEAL:: grads[q]= 2.00 4.00 2.00 4.00 -DEAL:: curls[q]= -2.00 -DEAL:: grads[q]= 2.00 4.00 2.00 4.00 -DEAL::FE=FESystem<3>[FE_Q<3>(1)^3] -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 -DEAL:: curls[q]= -6.00 9.00 -3.00 -DEAL:: grads[q]= 3.00 6.00 12.0 3.00 6.00 12.0 3.00 6.00 12.0 +DEAL::FE=FESystem<2>[FE_Q<2>(2)^2] +DEAL:: curls[q]= 0.423 +DEAL:: grads[q]= 0.00 0.00 0.423 8.33e-16 +DEAL:: curls[q]= 1.58 +DEAL:: grads[q]= 0.00 0.00 1.58 -1.11e-16 +DEAL:: curls[q]= 0.423 +DEAL:: grads[q]= 0.00 0.00 0.423 -6.66e-16 +DEAL:: curls[q]= 1.58 +DEAL:: grads[q]= 0.00 0.00 1.58 2.22e-16 -- 2.39.5