/**
* Compute the error for the solution of a system.
* See the other #integrate_difference#.
+ *
+ * The additional argument #weight# allows
+ * to evaluate weighted norms. This is useful
+ * for weighting the error of different parts
+ * differently. A special use is
+ * to have #weight=0# in some parts of the
+ * domain, e.g. at
+ * the location of a shock and #weight=1#
+ * elsewhere. This allows convergence tests
+ * in smooth parts of in general discontinuous
+ * solutions.
+ * By default, no weighting function is given,
+ * i.e. weight=1 in the whole domain.
*/
static void integrate_difference (const DoFHandler<dim> &dof,
const Vector<double> &fe_function,
const VectorFunction<dim>&exact_solution,
Vector<float> &difference,
const Quadrature<dim> &q,
- const NormType &norm);
+ const NormType &norm,
+ const Function<dim> *weight=0);
/**
const VectorFunction<dim>&exact_solution,
Vector<float> &difference,
const Quadrature<dim> &q,
- const NormType &norm)
+ const NormType &norm,
+ const Function<dim> *weight)
{
Assert(norm != mean , ExcNotUseful());
Assert (false, ExcNotImplemented());
};
+ // now weight the values
+ // at the quadrature points due
+ // to the weighting function
+ if (weight)
+ {
+ vector<double> w(n_q_points);
+ weight->value_list(fe_values.get_quadrature_points(),w);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ psi_scalar[q]*=w[q];
+ }
+
// ok, now we have the integrand,
// let's compute the integral,
// which is
- // sum_j psi_j JxW_j
+ // sum_j psi_j weight_j JxW_j
// (or |psi_j| or |psi_j|^2
switch (norm)
{
for (unsigned int k=0; k<fe.n_components; ++k)
psi_square[q] += sqr_point(psi[q][k]);
+ // now weight the values
+ // at the quadrature points due
+ // to the weighting function
+ if (weight)
+ {
+ vector<double> w(n_q_points);
+ weight->value_list(fe_values.get_quadrature_points(),w);
+ for (unsigned int q=0; q<n_q_points; ++q)
+ psi_square[q]*=w[q];
+ }
+
// add seminorm to L_2 norm or
// to zero
diff += inner_product (psi_square.begin(), psi_square.end(),