]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
prepared for lp-norms
authorguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 1 Aug 2002 08:20:10 +0000 (08:20 +0000)
committerguido <guido@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 1 Aug 2002 08:20:10 +0000 (08:20 +0000)
git-svn-id: https://svn.dealii.org/trunk@6308 0785d39b-7218-0410-832d-ea1e28bc413d

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

index 62b20af51ab1df7453098956930bb06c86ec3d90..160619dc763e30279cf5f3584a53d321d29fcae5 100644 (file)
@@ -305,10 +305,17 @@ class VectorTools
                                      *  absolute value of the
                                      *  function is integrated.
                                      * @item @p{L2_norm}: the square
-                                     *    of the function is
-                                     *    integrated on each cell;
-                                     *    afterwards the root is
-                                     *    taken of this value.
+                                     * of the function is integrated
+                                     * and the the square root of the
+                                     * result is computed on each
+                                     * cell.
+                                     * @item @p{Lp_norm}: the
+                                     * absolute value to the pth
+                                     * power is integrated and the
+                                     * pth root is computed on each
+                                     * cell. The exponent @p{p} is
+                                     * the last parameter of the
+                                     * function.
                                      *  @item @p{Linfty_norm}: the
                                      *  maximum absolute value of the
                                      *  function.
@@ -316,8 +323,11 @@ class VectorTools
                                      *    square of the function
                                      *    gradient is integrated on
                                      *    each cell; afterwards the
-                                     *    root is taken of this *
+                                     *    root is taken of this
                                      *    value.
+                                     * @item @p{W1p_seminorm}: this
+                                     * is the @p{Lp_norm} of the
+                                     * gradient.
                                      * @item @p{H1_norm}: the square
                                      *    of the function plus the
                                      *    square of the function
@@ -330,14 +340,21 @@ class VectorTools
                                      *    square of the
                                      *    @p{H1_seminorm}.
                                      *  @end{itemize}
-                                     */
+                                     * @item @p{W1p_norm}: like
+                                     * @p{H1_norm}, but for
+                                     * @p{Lp_norm} instead of
+                                     * @p{L2_norm}
+                                    */
     enum NormType {
          mean,
          L1_norm,
          L2_norm,
          Linfty_norm,
          H1_seminorm,
-         H1_norm
+         H1_norm,
+         Lp_norm,
+         W1p_seminorm,
+         W1p_norm
     };
     
                                     /**
@@ -736,18 +753,13 @@ class VectorTools
                                      * Compute the error of the finite element solution.
                                      * Integrate the difference between
                                      * a finite element function and
-                                     * the reference function, which
+                                     * a reference function, which
                                      * is given as a continuous function
                                      * object.
                                      *
-                                     * Note that this function returns
-                                     * its results in a vector of @p{float}s,
-                                     * rather than in a vector of @p{double}s,
-                                     * since accuracy is not that important
-                                     * here and to save memory. During
-                                     * computation of the results, the full
-                                     * accuracy of the @p{double} data type is
-                                     * used.
+                                     * The value of @p{exponent} is
+                                     * used for computing $L^p$-norms
+                                     * and $W^{1,p}$-norms.
                                      *
                                      * The additional argument
                                      * @p{weight} allows to evaluate
@@ -790,7 +802,8 @@ class VectorTools
                                      OutVector             &difference,
                                      const Quadrature<dim> &q,
                                      const NormType        &norm,
-                                     const Function<dim>   *weight=0);
+                                     const Function<dim>   *weight=0,
+                                     const double exponent = 2.);
 
                                     /**
                                      * Calls the @p{integrate_difference}
@@ -804,7 +817,8 @@ class VectorTools
                                      OutVector             &difference,
                                      const Quadrature<dim> &q,
                                      const NormType        &norm,
-                                     const Function<dim>   *weight=0);
+                                     const Function<dim>   *weight=0,
+                                     const double exponent = 2.);
 
                                     /**
                                      * Mean-value filter for Stokes.
index 53efe8b20d37ff6459360637a96cf4d393eee4df..821fbdc10aac05446248c619ecc7dc58a64b65c7 100644 (file)
@@ -1058,7 +1058,8 @@ VectorTools::integrate_difference (const Mapping<dim>    &mapping,
                                   OutVector             &difference,
                                   const Quadrature<dim> &q,
                                   const NormType        &norm,
-                                  const Function<dim>   *weight)
+                                  const Function<dim>   *weight,
+                                  double exponent)
 {
   const unsigned int        n_q_points   = q.n_quadrature_points;
   const FiniteElement<dim> &fe           = dof.get_fe();
@@ -1070,16 +1071,39 @@ VectorTools::integrate_difference (const Mapping<dim>    &mapping,
       Assert ((weight->n_components==1) || (weight->n_components==n_components),
              ExcDimensionMismatch(weight->n_components, n_components));
     }
-  
+
   difference.reinit (dof.get_tria().n_active_cells());
   
+  switch (norm)
+    {
+      case L2_norm:
+      case H1_seminorm:
+      case H1_norm:
+       exponent = 2.;
+       break;
+      case L1_norm:
+       exponent = 1.;
+       break;
+      default:
+       break;
+    }
+  
   UpdateFlags update_flags = UpdateFlags (update_q_points  |
                                          update_JxW_values);
-  if (norm != H1_seminorm)
-    update_flags = UpdateFlags(update_flags | update_values);
-  
-  if ((norm==H1_seminorm) || (norm==H1_norm))
-    update_flags = UpdateFlags (update_flags | update_gradients);
+  switch (norm)
+    {
+      case H1_seminorm:
+      case W1p_seminorm:
+       update_flags |= UpdateFlags (update_gradients);
+       break;
+      case H1_norm:
+      case W1p_norm:
+       update_flags |= UpdateFlags (update_gradients);
+                                        // no break!
+      default:
+       update_flags |= UpdateFlags (update_values);
+       break;
+    }  
   
   FEValues<dim> fe_values(mapping, fe, q, update_flags);
 
@@ -1115,267 +1139,247 @@ VectorTools::integrate_difference (const Mapping<dim>    &mapping,
        {
          if (weight->n_components>1)
            weight->vector_value_list (fe_values.get_quadrature_points(),
-                                     weight_vectors);
+                                      weight_vectors);
          else
-           weight->value_list (fe_values.get_quadrature_points(),
-                              weight_values);
+           {
+             weight->value_list (fe_values.get_quadrature_points(),
+                                 weight_values);
+//           for (unsigned int k=0;k<n_q_points;++k)
+//             weight_vectors[k] = weight_values[k];
+           }
+       }
+      
+      if (update_flags & update_values)
+       {
+                                          // first compute the exact solution
+                                          // (vectors) at the quadrature points
+                                          // try to do this as efficient as
+                                          // possible by avoiding a second
+                                          // virtual function call in case
+                                          // the function really has only
+                                          // one component
+         if (fe_is_system)
+           exact_solution.vector_value_list (fe_values.get_quadrature_points(),
+                                             psi_values);
+         else
+           {
+             exact_solution.value_list (fe_values.get_quadrature_points(),
+                                        tmp_values);
+             for (unsigned int i=0; i<n_q_points; ++i)
+               psi_values[i](0) = tmp_values[i];
+           }
+         
+                                          // then subtract finite element
+                                          // fe_function
+         fe_values.get_function_values (fe_function, function_values);
+         for (unsigned int q=0; q<n_q_points; ++q)
+           psi_values[q] -= function_values[q];
+       }
+
+                                      // Do the same for gradients, if required
+      if (update_flags & update_gradients)
+       {
+                                          // try to be a little clever
+                                          // to avoid recursive virtual
+                                          // function calls when calling
+                                          // @p{gradient_list} for functions
+                                          // that are really scalar
+                                          // functions
+         if (fe_is_system)
+           exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
+                                                psi_grads);
+         else
+           {
+             exact_solution.gradient_list (fe_values.get_quadrature_points(),
+                                           tmp_gradients);
+             for (unsigned int i=0; i<n_q_points; ++i)
+               psi_grads[i][0] = tmp_gradients[i];
+           }      
+         
+                                          // then subtract finite element
+                                          // function_grads
+         fe_values.get_function_grads (fe_function, function_grads);
+         for (unsigned int k=0; k<n_components; ++k)
+           for (unsigned int q=0; q<n_q_points; ++q)
+             psi_grads[q][k] -= function_grads[q][k];
        }
       
       switch (norm)
-       {
-         case mean:
-         case L1_norm:
-         case L2_norm:
-         case Linfty_norm:
-         case H1_norm:
-         {
-                                            // first compute the exact solution
-                                            // (vectors) at the quadrature points
-                                            // try to do this as efficient as
-                                            // possible by avoiding a second
-                                            // virtual function call in case
-                                            // the function really has only
-                                            // one component
-           if (fe_is_system)
-             exact_solution.vector_value_list (fe_values.get_quadrature_points(),
-                                               psi_values);
-           else
+       {
+         case mean:
+           if (weight != 0)
              {
-               exact_solution.value_list (fe_values.get_quadrature_points(),
-                                          tmp_values);
-               for (unsigned int i=0; i<n_q_points; ++i)
-                 psi_values[i](0) = tmp_values[i];
-             };
+                                                // Different weights for
+                                                // each component?
+               if (weight->n_components > 1)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                     psi_scalar[q] = 0;
+                                                      // weighted mean value
+                     for (unsigned int i=0; i<n_components; ++i)
+                       psi_scalar[q] += psi_values[q](i)
+                                        * weight_vectors[q](i);
+                   }
+               else // weight->n_components == 1
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                     psi_scalar[q] = 0;
+                                                      // weighted mean value
+                     for (unsigned int i=0; i<n_components; ++i)
+                       psi_scalar[q] += psi_values[q](i)
+                                        * weight_values[q];
+                   }
+             }
+           else // no weight function
+             for (unsigned int q=0; q<n_q_points; ++q)
+               {
+                 psi_scalar[q] = 0;
+                                                  // mean value
+                 for (unsigned int i=0; i<n_components; ++i)
+                   psi_scalar[q] += psi_values[q](i);
+               }
            
-                                            // then subtract finite element
-                                            // fe_function
-           fe_values.get_function_values (fe_function, function_values);
-           for (unsigned int q=0; q<n_q_points; ++q)
-             psi_values[q] -= function_values[q];
-
-           switch (norm)
+                                            // Integration on one cell
+           diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+                                      fe_values.get_JxW_values().begin(),
+                                      0.0);
+           break;
+                                            // Compute (weighted) squares
+                                            // in each quadrature point
+         case L2_norm:
+         case H1_norm:
+           if (weight != 0)
              {
-               case mean:
-                     if (weight != 0)
-                       {
-                                                          // Different weights for
-                                                          // each component?
-                         if (weight->n_components > 1)
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             {
-                               psi_scalar[q] = 0;
-                               // weighted mean value
-                               for (unsigned int i=0; i<n_components; ++i)
-                                 psi_scalar[q] += psi_values[q](i)
-                                                  * weight_vectors[q](i);
-                             }
-                         else // weight->n_components == 1
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             {
-                               psi_scalar[q] = 0;
-                               // weighted mean value
-                               for (unsigned int i=0; i<n_components; ++i)
-                                 psi_scalar[q] += psi_values[q](i)
-                                                  * weight_values[q];
-                             }
-                       }
-                     else // no weight function
-                       for (unsigned int q=0; q<n_q_points; ++q)
-                         {
-                           psi_scalar[q] = 0;
-                               // mean value
-                           for (unsigned int i=0; i<n_components; ++i)
-                             psi_scalar[q] += psi_values[q](i);
-                         }
-
-                                                      // Integration on one cell
-                     diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-                                                fe_values.get_JxW_values().begin(),
-                                                0.0);
-                     break;
-                                                      // Compute (weighted) squares
-                                                      // in each quadrature point
-               case L2_norm:
-               case H1_norm:
-                     if (weight != 0)
-                       {
-                                                          // Different weights for
-                                                          // each component?
-                         if (weight->n_components > 1)
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             {
-                               psi_scalar[q] = 0;
-                               // weighted scalar product
-                               for (unsigned int i=0; i<n_components; ++i)
-                                 psi_scalar[q] += psi_values[q](i)*psi_values[q](i)
-                                                  * weight_vectors[q](i);
-                             }
-                         else // weight->n_components == 1
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             psi_scalar[q] = psi_values[q].norm_sqr()
-                                             * weight_values[q];
-                       }
-                     else // no weight function
-                       for (unsigned int q=0; q<n_q_points; ++q)
-                         psi_scalar[q] = psi_values[q].norm_sqr();
-
-                                                      // Integration on one cell
-                     diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-                                                fe_values.get_JxW_values().begin(),
-                                                0.0);
-                     if (norm == L2_norm)
-                       diff=std::sqrt(diff);
-                     break;
-               case L1_norm:
-                     if (weight != 0)
-                       {
-                                                          // Different weights for
-                                                          // each component?
-                         if (weight->n_components > 1)
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             {
-                               psi_scalar[q] = 0;
-                               // weighted scalar product
-                               for (unsigned int i=0; i<n_components; ++i)
-                                 psi_scalar[q] += std::fabs(psi_values[q](i))
-                                                  * weight_vectors[q](i);
-                             }
-                         else // weight->n_components == 1
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             psi_scalar[q] = psi_values[q].l1_norm()
-                                             * weight_values[q];
-                       }
-                     else // no weight function
-                       for (unsigned int q=0; q<n_q_points; ++q)
-                         psi_scalar[q] = psi_values[q].l1_norm();                    
-
-                                                      // Integration on one cell
-                     diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-                                                fe_values.get_JxW_values().begin(),
-                                                0.0);
-                     if (norm == L2_norm)
-                       diff=std::sqrt(diff);
-                     break;
-               case Linfty_norm:
-                     if (weight != 0)
+                                                // Different weights for
+                                                // each component?
+               if (weight->n_components > 1)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                     psi_scalar[q] = 0;
+                                                      // weighted scalar product
+                     for (unsigned int i=0; i<n_components; ++i)
+                       psi_scalar[q] += psi_values[q](i)*psi_values[q](i)
+                                        * weight_vectors[q](i);
+                   }
+               else // weight->n_components == 1
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_scalar[q] = psi_values[q].norm_sqr()
+                                   * weight_values[q];
+             }
+           else // no weight function
+             for (unsigned int q=0; q<n_q_points; ++q)
+               psi_scalar[q] = psi_values[q].norm_sqr();
+           
+                                            // Integration on one cell
+           diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+                                      fe_values.get_JxW_values().begin(),
+                                      0.0);
+           if (norm == L2_norm)
+             diff=std::sqrt(diff);
+           break;
+         case L1_norm:
+           if (weight != 0)
+             {
+                                                // Different weights for
+                                                // each component?
+               if (weight->n_components > 1)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                     psi_scalar[q] = 0;
+                                                      // weighted scalar product
+                     for (unsigned int i=0; i<n_components; ++i)
+                       psi_scalar[q] += std::fabs(psi_values[q](i))
+                                        * weight_vectors[q](i);
+                   }
+               else // weight->n_components == 1
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_scalar[q] = psi_values[q].l1_norm()
+                                   * weight_values[q];
+             }
+           else // no weight function
+             for (unsigned int q=0; q<n_q_points; ++q)
+               psi_scalar[q] = psi_values[q].l1_norm();                      
+           
+                                            // Integration on one cell
+           diff = std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+                                      fe_values.get_JxW_values().begin(),
+                                      0.0);
+           if (norm == L2_norm)
+             diff=std::sqrt(diff);
+           break;
+         case Linfty_norm:
+           if (weight != 0)
+             {
+                                                // Different weights for
+                                                // each component?
+               if (weight->n_components > 1)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                     psi_scalar[q] = 0;
+                     for (unsigned int i=0; i<n_components; ++i)
                        {
-                                                          // Different weights for
-                                                          // each component?
-                         if (weight->n_components > 1)
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             {
-                               psi_scalar[q] = 0;
-                               for (unsigned int i=0; i<n_components; ++i)
-                                 {
-                                   double newval = std::fabs(psi_values[q](i))
-                                                   * weight_vectors[q](i);
-                                   if (psi_scalar[q]<newval)
-                                     psi_scalar[q] = newval;
-                                 }
-                             }
-                         else // weight->n_components == 1
-                           for (unsigned int q=0; q<n_q_points; ++q)
-                             psi_scalar[q] = psi_values[q].linfty_norm()
-                                             * weight_values[q];
+                         double newval = std::fabs(psi_values[q](i))
+                                         * weight_vectors[q](i);
+                         if (psi_scalar[q]<newval)
+                           psi_scalar[q] = newval;
                        }
-                     else // no weight function
-                       for (unsigned int q=0; q<n_q_points; ++q)
-                         psi_scalar[q] = psi_values[q].linfty_norm();
-
-                                                      // Maximum on one cell
-                     diff = *std::max_element (psi_scalar.begin(), psi_scalar.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!
-                                            // Until now, @p{diff} includes the
-                                            // square of the L2_norm.
-
-                                            // in praxi: first compute
-                                            // exact fe_function vector
-                                            //
-                                            // try to be a little clever
-                                            // to avoid recursive virtual
-                                            // function calls when calling
-                                            // @p{gradient_list} for functions
-                                            // that are really scalar
-                                            // functions
-           if (fe_is_system)
-             exact_solution.vector_gradient_list (fe_values.get_quadrature_points(),
-                                                  psi_grads);
-           else
+                   }
+               else // weight->n_components == 1
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_scalar[q] = psi_values[q].linfty_norm()
+                                   * weight_values[q];
+             }
+           else // no weight function
+             for (unsigned int q=0; q<n_q_points; ++q)
+               psi_scalar[q] = psi_values[q].linfty_norm();
+           
+                                            // Maximum on one cell
+           diff = *std::max_element (psi_scalar.begin(), psi_scalar.end());
+           break;
+         case H1_seminorm:
+           break;
+         default:
+           Assert (false, ExcNotImplemented());
+       }
+      
+      if (norm == H1_seminorm || norm == H1_norm)
+       {
+                                          // take square of integrand
+         std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
+         
+         for (unsigned int k=0; k<n_components; ++k)
+           if (weight != 0)
              {
-               exact_solution.gradient_list (fe_values.get_quadrature_points(),
-                                             tmp_gradients);
-               for (unsigned int i=0; i<n_q_points; ++i)
-                 psi_grads[i][0] = tmp_gradients[i];
-             };
-
-                                            // then subtract finite element
-                                            // function_grads
-           fe_values.get_function_grads (fe_function, function_grads);
-           for (unsigned int k=0; k<n_components; ++k)
+                                                // Different weights for
+                                                // each component?
+               if (weight->n_components > 1)
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   {
+                                                      // weighted scalar product
+                     psi_scalar[q] += sqr_point(psi_grads[q][k])
+                                      * weight_vectors[q](k);
+                   }
+               else // weight->n_components == 1
+                 for (unsigned int q=0; q<n_q_points; ++q)
+                   psi_scalar[q] = sqr_point(psi_grads[q][k])
+                                   * weight_values[q];
+             }
+           else // no weight function
              for (unsigned int q=0; q<n_q_points; ++q)
-               psi_grads[q][k] -= function_grads[q][k];
-
-                                            // take square of integrand
-           std::fill_n (psi_scalar.begin(), n_q_points, 0.0);
-
-           for (unsigned int k=0; k<n_components; ++k)
-             if (weight != 0)
-               {
-                                                  // Different weights for
-                                                  // each component?
-                 if (weight->n_components > 1)
-                   for (unsigned int q=0; q<n_q_points; ++q)
-                     {
-                               // weighted scalar product
-                       psi_scalar[q] += sqr_point(psi_grads[q][k])
-                                        * weight_vectors[q](k);
-                     }
-                 else // weight->n_components == 1
-                   for (unsigned int q=0; q<n_q_points; ++q)
-                     psi_scalar[q] = sqr_point(psi_grads[q][k])
-                                     * weight_values[q];
-               }
-             else // no weight function
-               for (unsigned int q=0; q<n_q_points; ++q)
-                 psi_scalar[q] += sqr_point(psi_grads[q][k]);
-
+               psi_scalar[q] += sqr_point(psi_grads[q][k]);
+         
                                             // add seminorm to L_2 norm or
                                             // to zero
-           diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
-                                       fe_values.get_JxW_values().begin(),
-                                       0.0);
-           diff = std::sqrt(diff);
-
-           break;
-         };
-                                            
-         default:
-               Assert (false, ExcNotImplemented());
-       };
-
-
+         diff += std::inner_product (psi_scalar.begin(), psi_scalar.end(),
+                                     fe_values.get_JxW_values().begin(),
+                                     0.0);
+         diff = std::sqrt(diff);         
+       }
                                       // append result of this cell
                                       // to the end of the vector
       difference(index) = diff;
-    };
-};
+    }
+}
 
 
 template <int dim, class InVector, class OutVector>
@@ -1386,12 +1390,13 @@ VectorTools::integrate_difference (const DoFHandler<dim>    &dof,
                                   OutVector                &difference,
                                   const Quadrature<dim>    &q,
                                   const NormType           &norm,
-                                  const Function<dim>      *weight)
+                                  const Function<dim>      *weight,
+                                  double exponent)
 {
   Assert (DEAL_II_COMPAT_MAPPING, ExcCompatibility("mapping"));
   static const MappingQ1<dim> mapping;
   integrate_difference(mapping, dof, fe_function, exact_solution,
-                      difference, q, norm, weight);
+                      difference, q, norm, weight, exponent);
 }
 
 

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.