Fix a bug reported by Christoph Heiniger (including a description of how it should...

git-svn-id: https://svn.dealii.org/trunk@32394 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/doc/news/changes.h b/deal.II/doc/news/changes.h
index 56b56bf4dc7778ab27d0a9777f26e3ab46132513..be608b1720cd9f904691ff279c450378c5eb04a7 100644 (file)
--- a/deal.II/doc/news/changes.h
+++ b/deal.II/doc/news/changes.h
@@ -124,6 +124,12 @@ inconvenience this causes.
 <h3>Specific improvements</h3>
 
 <ol>
+  <li>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.
+  <br>
+  (Christoph Heiniger, Wolfgang Bangerth, 2014/02/03)
+
   <li>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 aaf9fe0a0a746d355bc9e6670484848e0bc5a565..163f03d2b7f0350ad49cb53e16aead969bea3dd3 100644 (file)
--- 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 74b1736bcdaf3f66f1106f2247fadbb4a53e184f..284011d5405f9c10c875329df24d2bf46aa10587 100644 (file)
--- 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<dim> tr;
   GridGenerator::hyper_cube(tr);
 
-  FE_Nedelec<dim> fe(3);
+  FESystem<dim> fe(FE_Q<dim>(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 9341da79437d2d4f559f0cfda2e31f4803b8d2a3..c3c1e10f0bf840396274d416df2eb906b677a623 100644 (file)
--- 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