From: bangerth Date: Fri, 15 Dec 2006 23:11:42 +0000 (+0000) Subject: Implement a somewhat working smoothness estimator. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=e5518050a73fe01eb829c95cc16b8046c015a263;p=dealii-svn.git Implement a somewhat working smoothness estimator. git-svn-id: https://svn.dealii.org/trunk@14259 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 9c547e13d0..08ce436718 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -263,18 +263,33 @@ estimate_smoothness (Vector &smoothness_indicators) const { const unsigned int N = 5; - const unsigned n_fourier_modes = int_pow (N+1, dim); - std::vector > k_vectors (n_fourier_modes); - for (unsigned int k=0; k > k_vectors; + std::vector abs_k_square; + + switch (dim) + { + case 2: { - case 2: - k_vectors[k][0] = deal_II_numbers::PI * (k / (N+1)); - k_vectors[k][1] = deal_II_numbers::PI * (k % (N+1)); - break; - default: - Assert (false, ExcNotImplemented()); + for (unsigned int i=0; i(deal_II_numbers::PI * i, + deal_II_numbers::PI * j)); + abs_k_square.push_back (k_times_k); + } + } + + break; } + + default: + Assert (false, ExcNotImplemented()); + } + + const unsigned n_fourier_modes = k_vectors.size(); QGauss<1> base_quadrature (2); QIterated quadrature (base_quadrature, N); @@ -321,27 +336,38 @@ estimate_smoothness (Vector &smoothness_indicators) const dof_values(i); } + // 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.); + 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 (transformed_values[f])); - F[1] += std::log (std::abs (transformed_values[f])) * - std::log (k_abs); + 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];