From b3e56941491f5ac2cf6905bee92579800f7aef9a Mon Sep 17 00:00:00 2001 From: bangerth Date: Fri, 15 Dec 2006 03:53:15 +0000 Subject: [PATCH] Make the program do something useful, and do it without bombing out... git-svn-id: https://svn.dealii.org/trunk@14251 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-27/step-27.cc | 223 ++++++++++++++++++++++++++-- 1 file changed, 214 insertions(+), 9 deletions(-) diff --git a/deal.II/examples/step-27/step-27.cc b/deal.II/examples/step-27/step-27.cc index 89deda111f..9c547e13d0 100644 --- a/deal.II/examples/step-27/step-27.cc +++ b/deal.II/examples/step-27/step-27.cc @@ -48,6 +48,8 @@ #include #include +#include + // Finally, this is as in previous // programs: using namespace dealii; @@ -67,6 +69,7 @@ class LaplaceProblem void solve (); void create_coarse_grid (); void refine_grid (); + void estimate_smoothness (Vector &smoothness_indicators) const; void output_results (const unsigned int cycle) const; Triangulation triangulation; @@ -74,6 +77,7 @@ class LaplaceProblem hp::DoFHandler dof_handler; hp::FECollection fe_collection; hp::QCollection quadrature_collection; + hp::QCollection face_quadrature_collection; ConstraintMatrix hanging_node_constraints; @@ -85,14 +89,41 @@ class LaplaceProblem }; + +template +class RightHandSide : public Function +{ + public: + RightHandSide () : Function () {}; + + virtual double value (const Point &p, + const unsigned int component) const; +}; + + +template +double +RightHandSide::value (const Point &p, + const unsigned int /*component*/) const +{ + double product = 1; + for (unsigned int d=0; d LaplaceProblem::LaplaceProblem () : dof_handler (triangulation) { - for (unsigned int degree=1; degree<5; ++degree) + for (unsigned int degree=2; degree<7; ++degree) { fe_collection.push_back (FE_Q(degree)); quadrature_collection.push_back (QGauss(degree+2)); + face_quadrature_collection.push_back (QGauss(degree+2)); } } @@ -136,6 +167,8 @@ void LaplaceProblem::assemble_system () update_values | update_gradients | update_q_points | update_JxW_values); + const RightHandSide rhs_function; + typename hp::DoFHandler::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); @@ -151,8 +184,12 @@ void LaplaceProblem::assemble_system () cell_rhs = 0; hp_fe_values.reinit (cell); - + const FEValues &fe_values = hp_fe_values.get_present_fe_values (); + + std::vector rhs_values (fe_values.n_quadrature_points); + rhs_function.value_list (fe_values.get_quadrature_points(), + rhs_values); for (unsigned int q_point=0; q_point::assemble_system () fe_values.JxW(q_point)); cell_rhs(i) += (fe_values.shape_value(i,q_point) * - 1.0 * + rhs_values[q_point] * fe_values.JxW(q_point)); } @@ -207,24 +244,170 @@ void LaplaceProblem::solve () hanging_node_constraints.distribute (solution); } + +unsigned int +int_pow (const unsigned int x, + const unsigned int n) +{ + unsigned int p=1; + for (unsigned int i=0; i +void +LaplaceProblem:: +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 base_quadrature (2); + QIterated quadrature (base_quadrature, N); + + std::vector > > + fourier_transforms (fe_collection.size()); + for (unsigned int fe=0; fe sum = 0; + for (unsigned int q=0; q x_q = quadrature.point(q); + sum += std::exp(std::complex(0,1) * + (k_vectors[k] * x_q)) * + fe_collection[fe].shape_value(i,x_q) * + quadrature.weight(q); + } + fourier_transforms[fe](k,i) = sum; + } + } + + 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); + + 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); + } + + // now try to fit a curve C|k|^{-s} to + // the fourier transform 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 void LaplaceProblem::refine_grid () { Vector estimated_error_per_cell (triangulation.n_active_cells()); - KellyErrorEstimator::estimate (dof_handler, - QGauss(3), + face_quadrature_collection, typename FunctionMap::type(), solution, estimated_error_per_cell); + Vector smoothness_indicators (triangulation.n_active_cells()); + estimate_smoothness (smoothness_indicators); + GridRefinement::refine_and_coarsen_fixed_number (triangulation, estimated_error_per_cell, 0.3, 0.03); + float max_smoothness = 0, + min_smoothness = smoothness_indicators.linfty_norm(); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + if (cell->refine_flag_set()) + { + max_smoothness = std::max (max_smoothness, + smoothness_indicators(index)); + min_smoothness = std::min (min_smoothness, + smoothness_indicators(index)); + } + } + const float cutoff_smoothness = (max_smoothness + min_smoothness) / 2; + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + if (cell->refine_flag_set() + && + (smoothness_indicators(index) > cutoff_smoothness)) + { + cell->clear_refine_flag(); + cell->set_active_fe_index (std::min (cell->active_fe_index() + 1, + fe_collection.size() - 1)); + } + } + triangulation.execute_coarsening_and_refinement (); } + + template void LaplaceProblem::output_results (const unsigned int cycle) const { @@ -241,17 +424,36 @@ void LaplaceProblem::output_results (const unsigned int cycle) const } { + Vector smoothness_indicators (triangulation.n_active_cells()); + estimate_smoothness (smoothness_indicators); + + + Vector fe_indices (triangulation.n_active_cells()); + { + typename hp::DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (unsigned int index=0; cell!=endc; ++cell, ++index) + { + + fe_indices(index) = cell->active_fe_index(); +// smoothness_indicators(index) *= std::sqrt(cell->diameter()); + } + } + const std::string filename = "solution-" + Utilities::int_to_string (cycle, 2) + - ".gnuplot"; + ".vtk"; DataOut > data_out; 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 (fe_indices, "fe_index"); data_out.build_patches (); std::ofstream output (filename.c_str()); - data_out.write_gnuplot (output); + data_out.write_vtk (output); } } @@ -323,7 +525,7 @@ void LaplaceProblem<2>::create_coarse_grid () triangulation.create_triangulation (vertices, cells, SubCellData()); - triangulation.refine_global (1); + triangulation.refine_global (3); } @@ -331,7 +533,7 @@ void LaplaceProblem<2>::create_coarse_grid () template void LaplaceProblem::run () { - for (unsigned int cycle=0; cycle<5; ++cycle) + for (unsigned int cycle=0; cycle<8; ++cycle) { std::cout << "Cycle " << cycle << ':' << std::endl; @@ -350,6 +552,9 @@ void LaplaceProblem::run () std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs() << std::endl; + std::cout << " Number of constraints : " + << hanging_node_constraints.n_constraints() + << std::endl; assemble_system (); solve (); -- 2.39.5