From: wolf Date: Wed, 14 Apr 2004 21:20:39 +0000 (+0000) Subject: Commit another one. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=791291987531331200187ff8d8e471860f0aeeb4;p=dealii-svn.git Commit another one. git-svn-id: https://svn.dealii.org/trunk@9000 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-15/step-15.cc b/deal.II/examples/step-15/step-15.cc index 00b33b7b93..ea997d0946 100644 --- a/deal.II/examples/step-15/step-15.cc +++ b/deal.II/examples/step-15/step-15.cc @@ -373,19 +373,108 @@ void MinimizationProblem<1>::initialize () -template -void MinimizationProblem::refine_grid () +template <> +void MinimizationProblem<1>::refine_grid () { + const unsigned int dim = 1; + Vector estimated_error_per_cell (triangulation.n_active_cells()); - typename FunctionMap::type neumann_boundary; + QTrapez quadrature; + FEValues fe_values (fe, quadrature, + update_values | update_gradients | + update_second_derivatives | + update_q_points | update_JxW_values); + + FEValues neighbor_fe_values (fe, quadrature, + update_gradients); + + std::vector local_values (quadrature.n_quadrature_points); + std::vector > local_gradients (quadrature.n_quadrature_points); + std::vector > local_2nd_derivs (quadrature.n_quadrature_points); + + DoFHandler::active_cell_iterator + cell = dof_handler.begin_active (), + endc = dof_handler.end (); + for (unsigned int index = 0; cell!=endc; ++cell, ++index) + { + fe_values.reinit (cell); + fe_values.get_function_values (present_solution, local_values); + fe_values.get_function_grads (present_solution, local_gradients); + fe_values.get_function_2nd_derivatives (present_solution, local_2nd_derivs); + + double cell_residual_norm = 0; + for (unsigned int q=0; qdiameter() * cell->diameter(); + + const double u_left = local_values[0]; + const double u_right = local_values[1]; + + const double u_prime_left = local_gradients[0][0]; + const double u_prime_right = local_gradients[1][0]; - KellyErrorEstimator::estimate (dof_handler, - Quadrature(1), - neumann_boundary, - present_solution, - estimated_error_per_cell); + const double x_left = fe_values.quadrature_point(0)[0]; + const double x_right = fe_values.quadrature_point(1)[0]; + Assert (x_left == cell->vertex(0)[0], ExcInternalError()); + Assert (x_right == cell->vertex(1)[0], ExcInternalError()); + + if (cell->at_boundary(0) == false) + { + DoFHandler::cell_iterator left_neighbor = cell->neighbor(0); + while (left_neighbor->has_children()) + left_neighbor = left_neighbor->child(1); + + neighbor_fe_values.reinit (left_neighbor); + neighbor_fe_values.get_function_grads (present_solution, local_gradients); + + const double neighbor_u_prime_left = local_gradients[1][0]; + + const double left_jump = std::pow(x_left-std::pow(u_left,3), 2) * + (std::pow(neighbor_u_prime_left,5) - + std::pow(u_prime_left,5)); + estimated_error_per_cell(index) += left_jump * left_jump * + cell->diameter(); + } + + if (cell->at_boundary(1) == false) + { + DoFHandler::cell_iterator right_neighbor = cell->neighbor(1); + while (right_neighbor->has_children()) + right_neighbor = right_neighbor->child(0); + + neighbor_fe_values.reinit (right_neighbor); + neighbor_fe_values.get_function_grads (present_solution, local_gradients); + + const double neighbor_u_prime_right = local_gradients[0][0]; + + const double right_jump = std::pow(x_right-std::pow(u_right,3), 2) * + (std::pow(neighbor_u_prime_right,5) - + std::pow(u_prime_right,5)); + estimated_error_per_cell(index) += right_jump * right_jump * + cell->diameter(); + } + + } + GridRefinement::refine_and_coarsen_fixed_number (triangulation, estimated_error_per_cell, 0.3, 0.03);