fix incorrect Hdiv computation

|| u ||_Hdiv
= || div u ||_L2
= sqrt( \int_\Omega (div u)^2 )
= sqrt( \int_\Omega (sum d/dx_i u_i)^2 )
and not
= sqrt( \int_\Omega (sum (d/dx_i u_i)^2) )
which was the old implementation.

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index 8b6d4c272b82cf8ecb4985fc9ad5627aa76374f7..454c44881c010038033f63d12322ab73da4ac3d3 100644 (file)
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -6329,9 +6329,8 @@ namespace VectorTools
               // take the trace of the derivatives, square it, multiply it
               // with the weight function
               for (unsigned int k=0; k<dim; ++k)
-                sum += (data.psi_grads[q][k][k] * data.psi_grads[q][k][k]) *
-                       data.weight_vectors[q](k);
-              diff += sum * fe_values.JxW(q);
+                sum += data.psi_grads[q][k][k] * data.weight_vectors[q](k);
+              diff += sum * sum * fe_values.JxW(q);
             }
           diff = std::sqrt(diff);
           break;