From: Wolfgang Bangerth <bangerth@colostate.edu>
Date: Thu, 29 Mar 2018 22:24:35 +0000 (-0600)
Subject: Fix VectorTools::integrate_difference() for NaN components.
X-Git-Tag: v9.0.0-rc1~200^2~3
X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=c29dd40ae14080107fdfdfaeddeca8a37e366ba3;p=dealii.git

Fix VectorTools::integrate_difference() for NaN components.

The VectorTools::integrate_difference() function allows users
to provide a weight function that can also serve as a component mask
to select individual components of the solution vector for error
computation. For components not selected, such a mask would then
simply be zero.

In some cases, the solution vector contains NaN numbers, for example
when one uses the FE_FaceQ element for certain components of the
solution vector and uses a quadrature formula for error evaluation
that has quadrature points in the interior of the cell. For any
regular solution component for which the component mask has a zero
weight, the value of that component will be multiplied by zero and
consequently does not add anything to the error computation. However,
if the NaNs of a FE_FaceQ are multiplied with zero weights, the result
is still a NaN, and adding it to the values times weights of the other
components results in NaNs -- in effect rendering it impossible to get
any information out of the VectorTools::integrate_difference()
function if one of the finite elements involved is FE_FaceQ.

This is now fixed by simply skipping vector components for which the
weight vector is zero. This has the same result as before for all
normal situations, but also properly skips the NaN case outlined
above.
---

diff --git a/include/deal.II/numerics/vector_tools.templates.h b/include/deal.II/numerics/vector_tools.templates.h
index cab72cbc41..008e03b53e 100644
--- a/include/deal.II/numerics/vector_tools.templates.h
+++ b/include/deal.II/numerics/vector_tools.templates.h
@@ -7089,7 +7089,8 @@ namespace VectorTools
             {
               Number sum = 0;
               for (unsigned int k=0; k<n_components; ++k)
-                sum += data.psi_values[q](k) * data.weight_vectors[q](k);
+                if (data.weight_vectors[q](k) != 0)
+                  sum += data.psi_values[q](k) * data.weight_vectors[q](k);
               diff_mean += sum * fe_values.JxW(q);
             }
           break;
@@ -7102,9 +7103,10 @@ namespace VectorTools
             {
               double sum = 0;
               for (unsigned int k=0; k<n_components; ++k)
-                sum += std::pow(
-                         static_cast<double>(numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k))),
-                         exponent/2.) * data.weight_vectors[q](k);
+                if (data.weight_vectors[q](k) != 0)
+                  sum += std::pow(static_cast<double>(numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k))),
+                                  exponent/2.) *
+                         data.weight_vectors[q](k);
               diff += sum * fe_values.JxW(q);
             }
 
@@ -7120,8 +7122,9 @@ namespace VectorTools
             {
               double sum = 0;
               for (unsigned int k=0; k<n_components; ++k)
-                sum += numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k)) *
-                       data.weight_vectors[q](k);
+                if (data.weight_vectors[q](k) != 0)
+                  sum += numbers::NumberTraits<Number>::abs_square(data.psi_values[q](k)) *
+                         data.weight_vectors[q](k);
               diff += sum * fe_values.JxW(q);
             }
           // Compute the root only, if no derivative values are added later
@@ -7133,8 +7136,9 @@ namespace VectorTools
         case W1infty_norm:
           for (unsigned int q=0; q<n_q_points; ++q)
             for (unsigned int k=0; k<n_components; ++k)
-              diff = std::max (diff, double(std::abs(data.psi_values[q](k)*
-                                                     data.weight_vectors[q](k))));
+              if (data.weight_vectors[q](k) != 0)
+                diff = std::max (diff, double(std::abs(data.psi_values[q](k)*
+                                                       data.weight_vectors[q](k))));
           break;
 
         case H1_seminorm:
@@ -7158,9 +7162,10 @@ namespace VectorTools
             {
               double sum = 0;
               for (unsigned int k=0; k<n_components; ++k)
-                sum += std::pow(
-                         data.psi_grads[q][k].norm_square(),
-                         exponent/2.) * data.weight_vectors[q](k);
+                if (data.weight_vectors[q](k) != 0)
+                  sum += std::pow(
+                           data.psi_grads[q][k].norm_square(),
+                           exponent/2.) * data.weight_vectors[q](k);
               diff += sum * fe_values.JxW(q);
             }
           diff = std::pow(diff, 1./exponent);
@@ -7172,8 +7177,9 @@ namespace VectorTools
             {
               double sum = 0;
               for (unsigned int k=0; k<n_components; ++k)
-                sum += data.psi_grads[q][k].norm_square() *
-                       data.weight_vectors[q](k);
+                if (data.weight_vectors[q](k) != 0)
+                  sum += data.psi_grads[q][k].norm_square() *
+                         data.weight_vectors[q](k);
               diff += sum * fe_values.JxW(q);
             }
           diff = std::sqrt(diff);
@@ -7197,7 +7203,8 @@ namespace VectorTools
               Number sum = 0;
               // take the trace of the derivatives scaled by the weight and square it
               for (unsigned int k=idx; k<idx+dim; ++k)
-                sum += data.psi_grads[q][k][k-idx] * std::sqrt(data.weight_vectors[q](k));
+                if (data.weight_vectors[q](k) != 0)
+                  sum += data.psi_grads[q][k][k-idx] * std::sqrt(data.weight_vectors[q](k));
               diff += numbers::NumberTraits<Number>::abs_square(sum) * fe_values.JxW(q);
             }
           diff = std::sqrt(diff);
@@ -7209,10 +7216,11 @@ namespace VectorTools
           double t = 0;
           for (unsigned int q=0; q<n_q_points; ++q)
             for (unsigned int k=0; k<n_components; ++k)
-              for (unsigned int d=0; d<dim; ++d)
-                t = std::max(t,
-                             double(std::abs(data.psi_grads[q][k][d]) *
-                                    data.weight_vectors[q](k)));
+              if (data.weight_vectors[q](k) != 0)
+                for (unsigned int d=0; d<dim; ++d)
+                  t = std::max(t,
+                               double(std::abs(data.psi_grads[q][k][d]) *
+                                      data.weight_vectors[q](k)));
 
           // then add seminorm to norm if that had previously been computed
           diff += t;