]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Add computation of H1 norms and seminorms.
authorWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 8 Apr 1998 13:19:38 +0000 (13:19 +0000)
committerWolfgang Bangerth <bangerth@math.tamu.edu>
Wed, 8 Apr 1998 13:19:38 +0000 (13:19 +0000)
git-svn-id: https://svn.dealii.org/trunk@159 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/deal.II/include/numerics/base.h
deal.II/deal.II/source/numerics/base.cc

index 93d808433d0fcdffaf40531c9e20a0d3613489ee..7acd072b3dc880f6ed35881fc7d2160a0d0d64eb 100644 (file)
@@ -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 <int dim>
 class ProblemBase {
index 993c0e6156df11fdc5032c17d41b25f6d5ed73f3..b96bb8ba933abc6d7563319a917c62bd8c04d605 100644 (file)
 
 
 
-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();
+};
 
 
 
@@ -162,6 +166,8 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &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<dim> fe_values(fe, q, update_flags);
   
                                   // loop over all cells
@@ -188,6 +194,7 @@ void ProblemBase<dim>::integrate_difference (const Function<dim>      &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<dim>::integrate_difference (const Function<dim>      &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<dim>::integrate_difference (const Function<dim>      &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<dim>::integrate_difference (const Function<dim>      &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<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());
        };

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.