]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Move things a bit forward.
authorbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 14 Feb 2007 19:34:42 +0000 (19:34 +0000)
committerbangerth <bangerth@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 14 Feb 2007 19:34:42 +0000 (19:34 +0000)
git-svn-id: https://svn.dealii.org/trunk@14473 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-27/doc/intro.dox
deal.II/examples/step-27/step-27.cc

index f594bbbdeaa797874e9eeefd4c82151d200b6736..648e2fc21b61e889518846009c3566d854d02f8f 100644 (file)
@@ -114,8 +114,16 @@ So what do we have to do to estimate the local smoothness of $u(x)$ on
 a cell $K$? Clearly, the first step is to compute the Fourier series
 of our solution. Fourier series being infinite series, we simplify our
 task by only computing the first few terms of the series, such that
-$|\vec k|\le N$ with a cut-off $N$. Computing this series is not
-particularly hard: from the definition
+$|\vec k|\le N$ with a cut-off $N$. (Let us parenthetically remark
+that we want to choose $N$ large enough so that we capture at least
+the variation of those shape functions that vary the most. On the
+other hand, we should not choose $N$ too large: clearly, a finite
+element function, being a polynomial, is in $C^\infty$ on any given
+cell, so the coefficients will have to decay exponentially at one
+point; since we want to estimate the smoothness of the function this
+polynomial approximates, not of the polynomial itself, we need to
+choose a reasonable cutoff for $N$.) Either way, computing this series
+is not particularly hard: from the definition
 @f[
        \hat U_{\vec k}
        = \frac 1{(2\pi)^{d/2}} \int_{\hat K} e^{i\vec k \cdot \vec x} \hat u(\hat x) dx
@@ -172,29 +180,29 @@ problem
        \min_{\beta,\mu}
        Q(\beta,\mu) = 
        \frac 12 \sum_{\vec k, |\vec k|\le N}
-       \left( \ln |\hat U_{\vec k}| - \beta + \mu |\vec k|\right)^2,
+       \left( \ln |\hat U_{\vec k}| - \beta + \mu \ln |\vec k|\right)^2,
 @f]
 where $\beta=\ln \alpha$. This is now a problem for which the
 optimality conditions $\frac{\partial Q}{\partial\beta}=0,
 \frac{\partial Q}{\partial\mu}=0$, are linear in $\beta,\mu$. We can
 write these conditions as follows:
 @f[
-       \begin{array}{cc}
+       \left(\begin{array}{cc}
        \sum_{\vec k, |\vec k|\le N} 1 &
        \sum_{\vec k, |\vec k|\le N} \ln |\vec k| 
        \\
        \sum_{\vec k, |\vec k|\le N} \ln |\vec k| &
        \sum_{\vec k, |\vec k|\le N} (\ln |\vec k|)^2 
-       \end{array}
-       \begin{array}{c}
+       \end{array}\right)
+       \left(\begin{array}{c}
        \beta \\ -\mu
-       \end{array}
+       \end{array}\right)
        =
-       \begin{array}{c}
+       \left(\begin{array}{c}
        \sum_{\vec k, |\vec k|\le N} \ln |\hat U_{\vec k}|
        \\
        \sum_{\vec k, |\vec k|\le N} \ln |\hat U_{\vec k}| \ln |\vec k| 
-       \end{array}
+       \end{array}\right)
 @f]
 This linear system is readily inverted to yield
 @f[
@@ -217,11 +225,16 @@ and
                 \left(\sum_{\vec k, |\vec k|\le N} (\ln |\vec k|)^2\right)
                -\left(\sum_{\vec k, |\vec k|\le N} \ln |\vec k|\right)^2}
        \left[
-         -
          \left(\sum_{\vec k, |\vec k|\le N} \ln |\vec k|\right)
          \left(\sum_{\vec k, |\vec k|\le N} \ln |\hat U_{\vec k}|\right)
-         +
+         -
          \left(\sum_{\vec k, |\vec k|\le N} 1\right)
          \left(\sum_{\vec k, |\vec k|\le N} \ln |\hat U_{\vec k}| \ln |\vec k| \right)
        \right].
 @f]
+
+While we are not particularly interested in the actual value of
+$\beta$, the formula above gives us a mean to calculate the value of
+the exponent $\mu$ that we can then use to determine that $\hat u(\hat
+x)$ is in $H^s(\hat K)$ with $s=\mu-\frac d2$.
+
index 08ce436718955e633ab99eee0d8257e8fa68fd01..3f7d1cae25f8b3fd68531ccce3faa1768d9b4b96 100644 (file)
@@ -4,7 +4,7 @@
 /*    $Id$       */
 /*    Version: $Name$                                          */
 /*                                                                */
-/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006 by the deal.II authors */
+/*    Copyright (C) 2000, 2001, 2002, 2003, 2004, 2006, 2007 by the deal.II authors */
 /*                                                                */
 /*    This file is subject to QPL and may not be  distributed     */
 /*    without copyright and license information. Please refer     */
@@ -261,26 +261,25 @@ void
 LaplaceProblem<dim>::
 estimate_smoothness (Vector<float> &smoothness_indicators) const
 {
-  const unsigned int N = 5;
+  const unsigned int N = 7;
 
+                                  // form all the Fourier vectors
+                                  // that we want to
+                                  // consider. exclude k=0 to avoid
+                                  // problems with |k|^{-mu} and also
+                                  // logarithms of |k|
   std::vector<Tensor<1,dim> > k_vectors;
-  std::vector<unsigned int>   abs_k_square;
-  
   switch (dim)
     {
       case 2:
       {
        for (unsigned int i=0; i<N; ++i)
          for (unsigned int j=0; j<N; ++j)
-           {
-             const unsigned int k_times_k = i*i + j*j;
-             if (k_times_k < N*N)
-               {
-                 k_vectors.push_back (Point<2>(deal_II_numbers::PI * i,
-                                               deal_II_numbers::PI * j));
-                 abs_k_square.push_back (k_times_k);
-               }
-           }
+           if (!((i==0) && (j==0))
+               &&
+               (i*i + j*j < N*N))
+             k_vectors.push_back (Point<2>(deal_II_numbers::PI * i,
+                                           deal_II_numbers::PI * j));
            
        break;
       }
@@ -290,16 +289,26 @@ estimate_smoothness (Vector<float> &smoothness_indicators) const
     }
 
   const unsigned n_fourier_modes = k_vectors.size();
+  std::vector<double> ln_k (n_fourier_modes);
+  for (unsigned int i=0; i<n_fourier_modes; ++i)
+    ln_k[i] = std::log (k_vectors[i].norm());
   
+
+                                  // assemble the matrices that do
+                                  // the Fourier transforms for each
+                                  // of the finite elements we deal
+                                  // with. note that these matrices
+                                  // are complex-valued, so we can't
+                                  // use FullMatrix
   QGauss<1>      base_quadrature (2);
   QIterated<dim> quadrature (base_quadrature, N);
   
   std::vector<Table<2,std::complex<double> > >
-    fourier_transforms (fe_collection.size());
+    fourier_transform_matrices (fe_collection.size());
   for (unsigned int fe=0; fe<fe_collection.size(); ++fe)
     {
-      fourier_transforms[fe].reinit (n_fourier_modes,
-                                    fe_collection[fe].dofs_per_cell);
+      fourier_transform_matrices[fe].reinit (n_fourier_modes,
+                                            fe_collection[fe].dofs_per_cell);
 
       for (unsigned int k=0; k<n_fourier_modes; ++k)
        for (unsigned int i=0; i<fe_collection[fe].dofs_per_cell; ++i)
@@ -313,68 +322,61 @@ estimate_smoothness (Vector<float> &smoothness_indicators) const
                       fe_collection[fe].shape_value(i,x_q) *
                       quadrature.weight(q);
              }
-           fourier_transforms[fe](k,i) = sum;
+           fourier_transform_matrices[fe](k,i)
+             = sum / std::pow(2*deal_II_numbers::PI, 1.*dim/2);
          }
     }
 
+                                  // the next thing is to loop over
+                                  // all cells and do our work there,
+                                  // i.e. to locally do the Fourier
+                                  // transform and estimate the decay
+                                  // coefficient
+  std::vector<std::complex<double> > fourier_coefficients (n_fourier_modes);
+  Vector<double>                     local_dof_values;
+  
   typename hp::DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
-  std::vector<std::complex<double> > transformed_values (n_fourier_modes);
   for (unsigned int index=0; cell!=endc; ++cell, ++index)
     {
-      Vector<double> dof_values (cell->get_fe().dofs_per_cell);
-      cell->get_dof_values (solution, dof_values);
+      local_dof_values.reinit (cell->get_fe().dofs_per_cell);
+      cell->get_dof_values (solution, local_dof_values);
 
+                                      // first compute the Fourier
+                                      // transform of the local
+                                      // solution
+      std::fill (fourier_coefficients.begin(), fourier_coefficients.end(), 0);
+      for (unsigned int f=0; f<n_fourier_modes; ++f)
+       for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
+         fourier_coefficients[f] += 
+           fourier_transform_matrices[cell->active_fe_index()](f,i)
+           *
+           local_dof_values(i);
+
+                                      // now we have to calculate the
+                                      // various contributions to the
+                                      // formula for mu
+      double            sum_1 = 0,
+                    sum_ln_k = 0,
+             sum_ln_k_square = 0,
+             sum_ln_U        = 0,
+             sum_ln_U_ln_k   = 0;
       for (unsigned int f=0; f<n_fourier_modes; ++f)
        {
-         transformed_values[f] = 0;
-         for (unsigned int i=0; i<cell->get_fe().dofs_per_cell; ++i)
-           transformed_values[f] +=
-             fourier_transforms[cell->active_fe_index()](f,i)
-             *
-             dof_values(i);
+         sum_1 += 1;
+         sum_ln_k += ln_k[f];
+         sum_ln_k_square += ln_k[f]*ln_k[f];
+         sum_ln_U += std::log (std::abs (fourier_coefficients[f]));
+         sum_ln_U_ln_k += std::log (std::abs (fourier_coefficients[f])) * ln_k[f];
        }
 
-                                      // for each abs_k value we have, find
-                                      // the largest fourier coefficient
-      std::vector<double>
-       max_fourier_coefficient (*max_element(abs_k_square.begin(),
-                                             abs_k_square.end()) + 1,
-                                0.);
-      for (unsigned int f=0; f<n_fourier_modes; ++f)
-       max_fourier_coefficient[abs_k_square[f]]
-         = std::max (max_fourier_coefficient[abs_k_square[f]],
-                     std::abs (transformed_values[f]));
-      
-                                      // now try to fit a curve C|k|^{-s} to
-                                      // the maximal fourier coefficients of
-                                      // the solution on this cell
-      double A[2][2] = { { 0,0 }, { 0,0 }};
-      double F[2]    = { 0,0 };
-
-      for (unsigned int f=0; f<max_fourier_coefficient.size(); ++f)
-       if (max_fourier_coefficient[f] > 0)
-         {
-           const double k_abs = std::sqrt(1.*f);
-
-           if (k_abs == 0)
-             continue;
-         
-           A[0][0] += 1;
-           A[1][0] += std::log (k_abs);
-           A[1][1] += std::pow (std::log (k_abs), 2.);
-
-           F[0] += std::log (std::abs (max_fourier_coefficient[f]));
-           F[1] += std::log (std::abs (max_fourier_coefficient[f])) *
-                   std::log (k_abs);
-       }
-      A[0][1] = A[1][0];
-      
-      const double det = A[0][0] * A[1][1] - A[0][1] * A[0][1];
-      const double s = (A[0][0]*F[1] - A[1][1]*F[0]) / det;
+      const double mu
+       = (1./(sum_1*sum_ln_k_square - sum_ln_k*sum_ln_k)
+          *
+          (sum_ln_k*sum_ln_U - sum_1*sum_ln_U_ln_k));
       
-      smoothness_indicators(index) = s;
+      smoothness_indicators(index) = mu - 1.*dim/2;
     }
 }
 
@@ -453,6 +455,10 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
     Vector<float> smoothness_indicators (triangulation.n_active_cells());
     estimate_smoothness (smoothness_indicators);
 
+    Vector<double> smoothness_field (dof_handler.n_dofs());
+    DoFTools::distribute_cell_to_dof_vector (dof_handler,
+                                            smoothness_indicators,
+                                            smoothness_field);
 
     Vector<float> fe_indices (triangulation.n_active_cells());
     {
@@ -474,7 +480,8 @@ void LaplaceProblem<dim>::output_results (const unsigned int cycle) const
 
     data_out.attach_dof_handler (dof_handler);
     data_out.add_data_vector (solution, "solution");
-    data_out.add_data_vector (smoothness_indicators, "smoothness");
+    data_out.add_data_vector (smoothness_indicators, "smoothness1");
+    data_out.add_data_vector (smoothness_field, "smoothness2");
     data_out.add_data_vector (fe_indices, "fe_index");
     data_out.build_patches ();
   

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.