Deal with the potential. Read parameters from file.

git-svn-id: https://svn.dealii.org/trunk@18951 0785d39b-7218-0410-832d-ea1e28bc413d

diff --git a/deal.II/examples/step-36/step-36.cc b/deal.II/examples/step-36/step-36.cc
index 9045ff7b0167f5dd965559dcd351dc8b966db07e..7282b86472e20ab7817c45b747ee1bcbb9961ccf 100644 (file)
--- a/deal.II/examples/step-36/step-36.cc
+++ b/deal.II/examples/step-36/step-36.cc
@@ -53,7 +53,7 @@ template <int dim>
 class EigenvalueProblem 
 {
   public:
-    EigenvalueProblem ();
+    EigenvalueProblem (const std::string &prm_file);
     void run ();
     
   private:
@@ -79,7 +79,7 @@ class EigenvalueProblem
                                  // @sect4{EigenvalueProblem::EigenvalueProblem}
 
 template <int dim>
-EigenvalueProblem<dim>::EigenvalueProblem ()
+EigenvalueProblem<dim>::EigenvalueProblem (const std::string &prm_file)
                :
                 fe (1),
                dof_handler (triangulation)
@@ -91,6 +91,8 @@ EigenvalueProblem<dim>::EigenvalueProblem ()
   parameters.declare_entry ("Potential", "0",
                            Patterns::Anything(),
                            "A functional description of the potential.");
+
+  parameters.read_input (prm_file);
 }
 
 
@@ -123,7 +125,7 @@ void EigenvalueProblem<dim>::assemble_system ()
 
   FEValues<dim> fe_values (fe, quadrature_formula, 
                           update_values   | update_gradients |
-                           update_JxW_values);
+                           update_quadrature_points | update_JxW_values);
 
   const unsigned int dofs_per_cell = fe.dofs_per_cell;
   const unsigned int n_q_points    = quadrature_formula.size();
@@ -139,6 +141,7 @@ void EigenvalueProblem<dim>::assemble_system ()
                         "x,y,z"),
                        parameters.get ("Potential"),
                        typename FunctionParser<dim>::ConstMap());
+  std::vector<double> potential_values (n_q_points);
 
   ConstraintMatrix constraints;
   VectorTools::interpolate_boundary_values (dof_handler,
@@ -156,13 +159,20 @@ void EigenvalueProblem<dim>::assemble_system ()
       cell_stiffness_matrix = 0;
       cell_mass_matrix = 0;
 
+      potential.value_list (fe_values.get_quadrature_points(),
+                           potential_values);
+      
       for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
        for (unsigned int i=0; i<dofs_per_cell; ++i)
          for (unsigned int j=0; j<dofs_per_cell; ++j)
            {
              cell_stiffness_matrix(i,j)
-               += (fe_values.shape_grad (i, q_point) *
-                   fe_values.shape_grad (j, q_point) *
+               += ((fe_values.shape_grad (i, q_point) *
+                    fe_values.shape_grad (j, q_point)
+                    +
+                    potential_values[q_point] *
+                    fe_values.shape_value (i, q_point) *
+                    fe_values.shape_value (j, q_point)) *
                    fe_values.JxW (q_point));
              cell_mass_matrix(i,j)
                += (fe_values.shape_value (i, q_point) *
@@ -238,8 +248,8 @@ int main (int argc, char **argv)
       PetscInitialize(&argc,&argv,0,0);
       deallog.depth_console (0);
       
-      EigenvalueProblem<2> laplace_problem_2d;
-      laplace_problem_2d.run ();
+      EigenvalueProblem<2> problem ("step-36.prm");
+      problem.run ();
     }
   catch (std::exception &exc)
     {