mean,
L1_norm,
L2_norm,
- Linfty_norm
+ Linfty_norm,
+ H1_seminorm,
+ H1_norm
};
of freedom, which can then be attached to a #DataOut# object to be printed.
Presently, there is the possibility to compute the following values from the
- difference, on each cell: #mean#, #L1_norm#, #L2_norm#, #Linfty_norm#.
+ difference, on each cell: #mean#, #L1_norm#, #L2_norm#, #Linfty_norm#,
+ #H1_seminorm#.
For the mean difference value, the reference function minus the numerical
solution is computed, not the other way round.
intended as a rule of thumb, rather they are though to illustrate that the
use of the wrong quadrature formula may show a significantly wrong result
and care should be taken to chose the right formula.
+
+ The $H_1$ seminorm is the $L_2$ norm of the gradient of the difference. The
+ full $H_1$ norm is the sum of the seminorm and the $L_2$ norm.
To get the {\it global} L_1 error, you have to sum up the entries in
#difference#, e.g. using #dVector::l1_norm# function.
To get the global mean difference, simply sum up the elements as above.
To get the L_\infty norm, take the maximum of the vector elements, e.g.
using the #dVector::linfty_norm# function.
+
+ For the global $H_1$ norm and seminorm, the same rule applies as for the
+ $L_2$ norm: compute the $l_2$ norm of the cell error vector.
*/
template <int dim>
class ProblemBase {
-inline double sqr (double x) {
+inline double sqr (const double x) {
return x*x;
};
+template <int dim>
+inline double sqr_point (const Point<dim> &p) {
+ return p.square();
+};
update_flags.q_points = true;
update_flags.jacobians = true;
update_flags.JxW_values = true;
+ if ((norm==H1_seminorm) || (norm==H1_norm))
+ update_flags.gradients = true;
FEValues<dim> fe_values(fe, q, update_flags);
// loop over all cells
case L1_norm:
case L2_norm:
case Linfty_norm:
+ case H1_norm:
{
// we need the finite element
// function \psi at the different
psi.begin(), ptr_fun(fabs));
break;
case L2_norm:
+ case H1_norm:
transform (psi.begin(), psi.end(),
psi.begin(), ptr_fun(sqr));
break;
+ default:
+ Assert (false, ExcNotImplemented());
};
// ok, now we have the integrand,
0.0);
break;
case L2_norm:
+ case H1_norm:
diff = sqrt(inner_product (psi.begin(), psi.end(),
fe_values.get_JxW_values().begin(),
0.0));
case Linfty_norm:
diff = *max_element (psi.begin(), psi.end());
break;
+ default:
+ Assert (false, ExcNotImplemented());
};
+
+ // note: the H1_norm uses the result
+ // of the L2_norm and control goes
+ // over to the next case statement!
+ if (norm != H1_norm)
+ break;
+ };
+
+ case H1_seminorm:
+ {
+ // note: the computation of the
+ // H1_norm starts at the previous
+ // case statement, but continues
+ // here!
+
+ // for H1_norm: re-square L2_norm.
+ diff = sqr(diff);
+
+ // same procedure as above, but now
+ // psi is a vector of gradients
+ const unsigned int n_dofs = fe.total_dofs;
+ const unsigned int n_q_points = q.n_quadrature_points;
+ const vector<vector<Point<dim> > > & shape_grads = fe_values.get_shape_grads();
+ vector<double> dof_values;
+ cell->get_dof_values (solution, dof_values);
+ vector<Point<dim> > psi;
+
+ // in praxi: first compute
+ // exact solution vector
+ exact_solution.gradient_list (fe_values.get_quadrature_points(),
+ psi);
+
+ // then subtract finite element
+ // solution
+ for (unsigned int j=0; j<n_q_points; ++j)
+ for (unsigned int i=0; i<n_dofs; ++i)
+ psi[j] -= dof_values[i]*shape_grads[i][j];
+
+ // take square of integrand
+ vector<double> psi_square (psi.size(), 0.0);
+ for (unsigned int i=0; i<n_q_points; ++i)
+ psi_square[i] = sqr_point(psi[i]);
+
+ // add seminorm to L_2 norm or
+ // to zero
+ diff += inner_product (psi_square.begin(), psi_square.end(),
+ fe_values.get_JxW_values().begin(),
+ 0.0);
+ diff = sqrt(diff);
+
break;
};
-
+
default:
Assert (false, ExcNotImplemented());
};