From: Wolfgang Bangerth 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(numbers::NumberTraits::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(numbers::NumberTraits::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::abs_square(data.psi_values[q](k)) * - data.weight_vectors[q](k); + if (data.weight_vectors[q](k) != 0) + sum += numbers::NumberTraits::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::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