From: Wolfgang Bangerth Date: Wed, 8 Apr 1998 13:19:38 +0000 (+0000) Subject: Add computation of H1 norms and seminorms. X-Git-Tag: v8.0.0~23095 X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=70da0f643c7fefc52207c6e2184fea4deca7dfa8;p=dealii.git Add computation of H1 norms and seminorms. git-svn-id: https://svn.dealii.org/trunk@159 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/deal.II/include/numerics/base.h b/deal.II/deal.II/include/numerics/base.h index 93d808433d..7acd072b3d 100644 --- a/deal.II/deal.II/include/numerics/base.h +++ b/deal.II/deal.II/include/numerics/base.h @@ -41,7 +41,9 @@ enum NormType { mean, L1_norm, L2_norm, - Linfty_norm + Linfty_norm, + H1_seminorm, + H1_norm }; @@ -176,7 +178,8 @@ enum NormType { 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. @@ -198,6 +201,9 @@ enum NormType { 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. @@ -210,6 +216,9 @@ enum NormType { 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 class ProblemBase { diff --git a/deal.II/deal.II/source/numerics/base.cc b/deal.II/deal.II/source/numerics/base.cc index 993c0e6156..b96bb8ba93 100644 --- a/deal.II/deal.II/source/numerics/base.cc +++ b/deal.II/deal.II/source/numerics/base.cc @@ -18,11 +18,15 @@ -inline double sqr (double x) { +inline double sqr (const double x) { return x*x; }; +template +inline double sqr_point (const Point &p) { + return p.square(); +}; @@ -162,6 +166,8 @@ void ProblemBase::integrate_difference (const Function &exact_sol 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 fe_values(fe, q, update_flags); // loop over all cells @@ -188,6 +194,7 @@ void ProblemBase::integrate_difference (const Function &exact_sol case L1_norm: case L2_norm: case Linfty_norm: + case H1_norm: { // we need the finite element // function \psi at the different @@ -240,9 +247,12 @@ void ProblemBase::integrate_difference (const Function &exact_sol 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, @@ -259,6 +269,7 @@ void ProblemBase::integrate_difference (const Function &exact_sol 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)); @@ -266,11 +277,63 @@ void ProblemBase::integrate_difference (const Function &exact_sol 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 > > & shape_grads = fe_values.get_shape_grads(); + vector dof_values; + cell->get_dof_values (solution, dof_values); + vector > 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 psi_square (psi.size(), 0.0); + for (unsigned int i=0; i