From 6b70e43ac8a309eb6c95580248d867ecb3e98421 Mon Sep 17 00:00:00 2001 From: bangerth Date: Wed, 14 Feb 2007 19:34:42 +0000 Subject: [PATCH] Move things a bit forward. git-svn-id: https://svn.dealii.org/trunk@14473 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-27/doc/intro.dox | 35 +++++-- deal.II/examples/step-27/step-27.cc | 137 +++++++++++++------------ 2 files changed, 96 insertions(+), 76 deletions(-) diff --git a/deal.II/examples/step-27/doc/intro.dox b/deal.II/examples/step-27/doc/intro.dox index f594bbbdea..648e2fc21b 100644 --- a/deal.II/examples/step-27/doc/intro.dox +++ b/deal.II/examples/step-27/doc/intro.dox @@ -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$. + diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 08ce436718..3f7d1cae25 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -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:: estimate_smoothness (Vector &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 > k_vectors; - std::vector abs_k_square; - switch (dim) { case 2: { for (unsigned int i=0; i(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 &smoothness_indicators) const } const unsigned n_fourier_modes = k_vectors.size(); + std::vector ln_k (n_fourier_modes); + for (unsigned int i=0; i base_quadrature (2); QIterated quadrature (base_quadrature, N); std::vector > > - fourier_transforms (fe_collection.size()); + fourier_transform_matrices (fe_collection.size()); for (unsigned int fe=0; fe &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 > fourier_coefficients (n_fourier_modes); + Vector local_dof_values; + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); - std::vector > transformed_values (n_fourier_modes); for (unsigned int index=0; cell!=endc; ++cell, ++index) { - Vector 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; fget_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; fget_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 - max_fourier_coefficient (*max_element(abs_k_square.begin(), - abs_k_square.end()) + 1, - 0.); - for (unsigned int f=0; 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::output_results (const unsigned int cycle) const Vector smoothness_indicators (triangulation.n_active_cells()); estimate_smoothness (smoothness_indicators); + Vector smoothness_field (dof_handler.n_dofs()); + DoFTools::distribute_cell_to_dof_vector (dof_handler, + smoothness_indicators, + smoothness_field); Vector fe_indices (triangulation.n_active_cells()); { @@ -474,7 +480,8 @@ void LaplaceProblem::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 (); -- 2.39.5