From 791291987531331200187ff8d8e471860f0aeeb4 Mon Sep 17 00:00:00 2001
From: wolf <wolf@0785d39b-7218-0410-832d-ea1e28bc413d>
Date: Wed, 14 Apr 2004 21:20:39 +0000
Subject: [PATCH] Commit another one.

git-svn-id: https://svn.dealii.org/trunk@9000 0785d39b-7218-0410-832d-ea1e28bc413d
---
 deal.II/examples/step-15/step-15.cc | 105 +++++++++++++++++++++++++---
 1 file changed, 97 insertions(+), 8 deletions(-)

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 <int dim>
-void MinimizationProblem<dim>::refine_grid ()
+template <>
+void MinimizationProblem<1>::refine_grid ()
 {
+  const unsigned int dim = 1;
+  
   Vector<float> estimated_error_per_cell (triangulation.n_active_cells());
 
-  typename FunctionMap<dim>::type neumann_boundary;
+  QTrapez<dim> quadrature;
+  FEValues<dim> fe_values (fe, quadrature,
+                           update_values   | update_gradients |
+                           update_second_derivatives |
+                           update_q_points | update_JxW_values);
+
+  FEValues<dim> neighbor_fe_values (fe, quadrature,
+                                    update_gradients);
+
+  std::vector<double> local_values (quadrature.n_quadrature_points);
+  std::vector<Tensor<1,dim> > local_gradients (quadrature.n_quadrature_points);
+  std::vector<Tensor<2,dim> > local_2nd_derivs (quadrature.n_quadrature_points);
+
+  DoFHandler<dim>::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; q<quadrature.n_quadrature_points; ++q)
+        {
+          const double x             = fe_values.quadrature_point(q)[0];
+          const double u             = local_values[q];
+          const double u_prime       = local_gradients[q][0];
+          const double u_doubleprime = local_2nd_derivs[q][0][0];
+          const double local_residual_value
+            = ((x-u*u*u) * std::pow(u_prime, 4) *
+               (u*u*u_prime*u_prime
+                +
+                5*(x-u*u*u)*u_doubleprime
+                +
+                2*u_prime*(1-3*u*u*u_prime)));
+          
+          cell_residual_norm += (local_residual_value * local_residual_value *
+                                 fe_values.JxW(q));
+        }
+
+      estimated_error_per_cell(index) = cell_residual_norm *
+                                        cell->diameter() * 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<dim>::estimate (dof_handler,
-				      Quadrature<dim-1>(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<dim>::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<dim>::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);
-- 
2.39.5