]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated step-29, now using new DataPostprocessor class to output data
authorallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 12 Sep 2007 23:39:11 +0000 (23:39 +0000)
committerallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Wed, 12 Sep 2007 23:39:11 +0000 (23:39 +0000)
git-svn-id: https://svn.dealii.org/trunk@15201 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-29/step-29.cc
deal.II/examples/step-29/step-29.prm

index cafd6323280546b52f4c7f71a491a785b6422811..2e6542348b4faf99e4d8863965d834b7c078d49a 100644 (file)
@@ -1,75 +1,81 @@
-/*    $Id: step-28.cc 14713 2007-05-27 04:05:41Z bangerth $       */
-/*    Version: $Name:  $                                          */
-/*                                                                */
-/*    Copyright (C) 2007 by the deal.II authors and Moritz Allmaras     */
-/*                                                                */
-/*    This file is subject to QPL and may not be  distributed     */
-/*    without copyright and license information. Please refer     */
-/*    to the file deal.II/doc/license.html for the  text  and     */
-/*    further information on this license.                        */
+                               // @sect3{Include files}
+
+                               // The following header files are unchanged 
+                               // from step-7 and have been discussed before:
 
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 #include <base/logstream.h>
-#include <base/parameter_handler.h>
-#include <base/subscriptor.h>
+
 #include <lac/vector.h>
 #include <lac/full_matrix.h>
 #include <lac/sparse_matrix.h>
-#include <lac/sparse_direct.h>
+
 #include <grid/tria.h>
 #include <grid/grid_generator.h>
 #include <grid/tria_accessor.h>
 #include <grid/tria_iterator.h>
 #include <grid/tria_boundary_lib.h>
+
 #include <dofs/dof_handler.h>
 #include <dofs/dof_accessor.h>
 #include <dofs/dof_tools.h>
-#include <fe/fe_system.h>
+
 #include <fe/fe_q.h>
 #include <fe/fe_values.h>
+
 #include <numerics/matrices.h>
 #include <numerics/data_out.h>
 #include <numerics/vectors.h>
+
 #include <fstream>
 
-#define DIM 2
+                               // This header file is needed for  
+                               // the ParameterHandler class that we will
+                               // use to read parameters from a 
+                               // configuration file during runtime:
+#include <base/parameter_handler.h>
 
+#include <lac/sparse_direct.h>
+#include <fe/fe_system.h>
+
+#define DIM 2
 
 using namespace dealii;
 
+
 template <int dim>
 class DirichletBoundaryValues : public Function<dim>
 {
   public:
     DirichletBoundaryValues() : Function<dim> (2) {};
-    
-    virtual void vector_value (        const Point<dim> &p,
-                                Vector<double>   &values) const;
-    virtual void vector_value_list (const std::vector<Point<dim> > &   points,
-                                    std::vector<Vector<double> > &             value_list) const;
+
+    virtual void vector_value (const Point<dim> &p,
+                              Vector<double>   &values) const;
+
+    virtual void vector_value_list (const std::vector<Point<dim> > &points,
+                                   std::vector<Vector<double> >   &value_list) const;
 };
 
 
 template <int dim>
 inline
-void DirichletBoundaryValues<dim>::vector_value (      const Point<dim> &      /*p*/,
-                                                       Vector<double> &                values) const 
+void DirichletBoundaryValues<dim>::vector_value (const Point<dim> &/*p*/,
+                                                Vector<double>   &values) const 
 {
   Assert (values.size() == 2, ExcDimensionMismatch (values.size(), 2));
-       
+
   values(0) = 1;
   values(1) = 0;
 }
 
 
 template <int dim>
-void DirichletBoundaryValues<dim>::vector_value_list (const std::vector<Point<dim> > & points,
-                                                     std::vector<Vector<double> > &            value_list) const 
+void DirichletBoundaryValues<dim>::vector_value_list (const std::vector<Point<dim> > &points,
+                                                     std::vector<Vector<double> >   &value_list) const 
 {
   Assert (value_list.size() == points.size(), 
-          ExcDimensionMismatch (value_list.size(), points.size()));
+         ExcDimensionMismatch (value_list.size(), points.size()));
 
   for (unsigned int p=0; p<points.size(); ++p)
     DirichletBoundaryValues<dim>::vector_value (points[p], value_list[p]);
@@ -81,53 +87,142 @@ class ParameterReader : public Subscriptor
   public:
     ParameterReader(ParameterHandler &);
     void read_parameters();
+
   private:
     void declare_parameters();
-    ParameterHandler *prm;
+    ParameterHandler &prm;
 };
 
 
 ParameterReader::ParameterReader(ParameterHandler &paramhandler)
-               :
-               prm(&paramhandler)
+  :
+  prm(paramhandler)
 {}
 
 
+void ParameterReader::read_parameters()
+{
+  declare_parameters();
+
+  const std::string parameter_file = "step-29.prm";
+  prm.read_input (parameter_file);
+}
+
+
 void ParameterReader::declare_parameters()
 {
-  prm->declare_entry(  "Number of refinements", "5",
-                       Patterns::Integer(1,10),
-                       "Number of global mesh refinement steps "
-                       "applied to initial coarse grid");
+  prm.enter_subsection ("Mesh & geometry parameters");
+
+    prm.declare_entry("Number of refinements", "6",
+                     Patterns::Integer(1,10),
+                     "Number of global mesh refinement steps "
+                     "applied to initial coarse grid");
+
+     prm.declare_entry("Focal distance", "0.3",
+                      Patterns::Double(0),
+                      "Distance of the focal point of the lens "
+                      "to the x-axis (or xy-plane in 3D)");
 
-  prm->declare_entry(  "Focal distance", "0.3",
-                       Patterns::Double(0),
-                       "Distance of the focal point of the lens "
-                       "to the x-axis (or xy-plane in 3D)");
+  prm.leave_subsection ();
 
-  prm->declare_entry(  "c", "1.5e5",
-                       Patterns::Double(),
-                       "Wave speed");
+  prm.enter_subsection ("Physical constants");
 
-  prm->declare_entry(  "omega", "1.5e7",
-                       Patterns::Double(),
-                       "Frequency");
+    prm.declare_entry("c", "1.5e5",
+                     Patterns::Double(),
+                     "Wave speed");
 
-  prm->declare_entry(  "Output file", "solution",
-                       Patterns::Anything(),
-                       "Name of the output file (without extension)");
+    prm.declare_entry("omega", "5.0e7",
+                     Patterns::Double(),
+                     "Frequency");
 
-  DataOutInterface<1>::declare_parameters (*prm);
+  prm.leave_subsection ();
+
+  prm.enter_subsection ("Output parameters");
+
+    prm.declare_entry("Output file", "solution",
+                     Patterns::Anything(),
+                     "Name of the output file (without extension)");
+
+    DataOutInterface<1>::declare_parameters (prm);
+
+  prm.leave_subsection ();
 }
 
 
-void ParameterReader::read_parameters()
+template <int dim>
+class Postprocessor : public DataPostprocessor<dim>
 {
-  declare_parameters();
+  public:
 
-  const std::string parameter_file = "step-29.prm";
-  prm->read_input (parameter_file);
+    using DataPostprocessor<dim>::compute_derived_quantities;
+
+    void compute_derived_quantities (
+                       std::vector< Vector< double > > &,
+                       const std::vector< Vector< double > > &, 
+                       const std::vector< std::vector< Tensor< 1, dim > > > &, 
+                       const std::vector< std::vector< Tensor< 2, dim > > > &, 
+                       const std::vector< Point< dim > >                    &
+                       ) const;
+
+    std::vector<std::string> get_names () const;
+    UpdateFlags              get_needed_update_flags () const;
+    unsigned int             n_output_variables () const;
+};
+
+
+template <int dim>
+std::vector<std::string>
+Postprocessor<dim>::get_names() const
+{
+  std::vector<std::string> field_names;
+
+  field_names.push_back("Re_u");
+  field_names.push_back("Im_u");
+  field_names.push_back("Intensity");
+
+  return field_names;
+}
+
+
+template <int dim>
+UpdateFlags
+Postprocessor<dim>::get_needed_update_flags () const
+{
+  return update_values;
+}
+
+
+template <int dim>
+unsigned int
+Postprocessor<dim>::n_output_variables () const
+{
+  return 3;
+}
+
+
+template <int dim>
+void
+Postprocessor<dim>::compute_derived_quantities (
+                       std::vector< Vector< double > >                       &computed_quantities,
+                       const std::vector< Vector< double > >                 &uh,
+                       const std::vector< std::vector< Tensor< 1, dim > > >  &/*duh*/,
+                       const std::vector< std::vector< Tensor< 2, dim > > >  &/*dduh*/,
+                       const std::vector< Point< dim > >                     &/*normals*/
+                       ) const
+{
+  Assert(computed_quantities.size() == uh.size(), 
+         ExcDimensionMismatch (computed_quantities.size(), uh.size()));
+
+  for (unsigned int i=0; i<computed_quantities.size(); i++)
+  {
+    Assert(computed_quantities[i].size() == 3, 
+           ExcDimensionMismatch (computed_quantities[i].size(), 3));
+    Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));
+
+    computed_quantities[i](0) = uh[i](0);
+    computed_quantities[i](1) = uh[i](1);
+    computed_quantities[i](2) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
+  }
 }
 
 
@@ -138,38 +233,63 @@ class UltrasoundProblem
     UltrasoundProblem (ParameterHandler &);
     ~UltrasoundProblem ();
     void run ();
-    
+
   private:
+    static double get_omega (ParameterHandler &);
+    static double get_c (ParameterHandler &);
+
     void make_grid ();
     void setup_system ();
     void assemble_system ();
     void solve ();
-    void postprocess ();
     void output_results () const;
 
-    ParameterHandler &prm; 
+    ParameterHandler       &prm; 
+
+    Triangulation<dim>     triangulation;
+    DoFHandler<dim>        dof_handler;
+    FESystem<dim>          fe;
 
-    Triangulation<dim>                 triangulation;
-    DoFHandler<dim>         dof_handler;
-    FESystem<dim>                                              fe;
-    
-    SparsityPattern         sparsity_pattern;
-    SparseMatrix<double>    system_matrix;
-               
-    Vector<double>          solution, system_rhs, intensity;
+    SparsityPattern        sparsity_pattern;
+    SparseMatrix<double>   system_matrix;      
+    Vector<double>         solution, system_rhs;
 
-    const double                                               c, omega;
+    const double           c, omega;
 };
 
 
 template <int dim>
-UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler&   param) 
-               :
-               prm(param),
-               dof_handler(triangulation),
-               fe(FE_Q<dim>(1), 2),
-               c(prm.get_double("c")), 
-               omega(prm.get_double("omega"))
+double
+UltrasoundProblem<dim>::get_omega (ParameterHandler &prm)
+{
+  prm.enter_subsection ("Physical constants");
+    double omega_tmp = prm.get_double("omega");
+  prm.leave_subsection ();
+
+  return omega_tmp;
+}
+
+
+template <int dim>
+double
+UltrasoundProblem<dim>::get_c (ParameterHandler &prm)
+{
+  prm.enter_subsection ("Physical constants");
+    double c_tmp = prm.get_double("c");
+  prm.leave_subsection ();
+
+  return c_tmp;
+}
+
+
+template <int dim>
+UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler&  param) 
+  :
+  prm(param),
+  dof_handler(triangulation),
+  fe(FE_Q<dim>(1), 2),
+  c(get_c(prm)),
+  omega(get_omega(prm))
 {}
 
 
@@ -183,33 +303,41 @@ UltrasoundProblem<dim>::~UltrasoundProblem ()
 template <int dim>
 void UltrasoundProblem<dim>::make_grid ()
 {
+  prm.enter_subsection ("Mesh & geometry parameters");
+
+  const double         focal_distance = prm.get_double("Focal distance");
+  const unsigned int   N_ref          = prm.get_integer("Number of refinements");
+
+  prm.leave_subsection ();
+
   GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
-       
+
   const Point<dim>     transducer = (dim == 2) ? 
                                     Point<dim> (0.5, 0.0) :
                                     Point<dim> (0.5, 0.5, 0.0), 
-                      focal_point = (dim == 2) ?
-                                    Point<dim> (0.5, prm.get_double("Focal distance")) :
-                                    Point<dim> (0.5, 0.5, prm.get_double("Focal distance"));
+                       focal_point = (dim == 2) ?
+                                     Point<dim> (0.5, focal_distance) :
+                                     Point<dim> (0.5, 0.5, focal_distance);
 
-  double radius = sqrt(        (focal_point.distance(transducer) * 
+  double radius = sqrt( (focal_point.distance(transducer) * 
                         focal_point.distance(transducer)) + 
                        ((dim==2) ? 0.01 : 0.02));
 
   typename Triangulation<dim>::cell_iterator
     cell = triangulation.begin (),
     endc = triangulation.end();
-         
+
   for (; cell!=endc; ++cell)
     for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
       if ( cell->face(face)->at_boundary() &&
           ((cell->face(face)->center() - transducer).square() < 0.01) )
-       cell->face(face)->set_boundary_indicator (1);
+
+        cell->face(face)->set_boundary_indicator (1);
 
   const HyperBallBoundary<dim> boundary(focal_point, radius);
   triangulation.set_boundary(1, boundary);
 
-  triangulation.refine_global (prm.get_integer("Number of refinements"));
+  triangulation.refine_global (N_ref);
 
   deallog << " Number of active cells:  "
          << triangulation.n_active_cells()
@@ -225,40 +353,39 @@ void UltrasoundProblem<dim>::setup_system ()
   dof_handler.distribute_dofs (fe);
 
   deallog << " Number of degrees of freedom: "
-          << dof_handler.n_dofs()
-          << std::endl;
+         << dof_handler.n_dofs()
+         << std::endl;
+
+  sparsity_pattern.reinit (dof_handler.n_dofs(),
+                          dof_handler.n_dofs(),
+                          dof_handler.max_couplings_between_dofs());
 
-  sparsity_pattern.reinit (    dof_handler.n_dofs(),
-                               dof_handler.n_dofs(),
-                               dof_handler.max_couplings_between_dofs());
   DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);
   sparsity_pattern.compress();
 
   system_matrix.reinit (sparsity_pattern);
   system_rhs.reinit (dof_handler.n_dofs());
   solution.reinit (dof_handler.n_dofs());
-  intensity.reinit(triangulation.n_active_cells());
 }
 
 
 template <int dim>
 void UltrasoundProblem<dim>::assemble_system () 
-{ 
+{
   const double om2 = omega * omega;
   const double c2 = c * c;
-       
-  QGauss<dim>          quadrature_formula(3);
-  QGauss<dim-1> face_quadrature_formula(3);
 
-  const unsigned int   n_q_points              = quadrature_formula.n_quadrature_points;
-  const unsigned int   n_face_q_points = face_quadrature_formula.n_quadrature_points;
+  QGauss<dim>    quadrature_formula(2);
+  QGauss<dim-1>  face_quadrature_formula(2);
 
-  const unsigned int   dofs_per_cell = fe.dofs_per_cell;
+  const unsigned int n_q_points              = quadrature_formula.n_quadrature_points,
+                    n_face_q_points  = face_quadrature_formula.n_quadrature_points,
+                    dofs_per_cell    = fe.dofs_per_cell;
 
-  FullMatrix<double>   cell_matrix (dofs_per_cell, dofs_per_cell);
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
 
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
-  
+
   FEValues<dim>  fe_values (fe, quadrature_formula, 
                            update_values | update_gradients |
                            update_JxW_values);
@@ -269,67 +396,70 @@ void UltrasoundProblem<dim>::assemble_system ()
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     endc = dof_handler.end();
+
   for (; cell!=endc; ++cell)
+  {
+    cell_matrix = 0;
+    fe_values.reinit (cell);
+
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
     {
-      cell_matrix = 0;
-      fe_values.reinit (cell);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       {
-         for (unsigned int j=0; j<dofs_per_cell; ++j)
-           {
-             if (fe.system_to_component_index(i).first == 
-                 fe.system_to_component_index(j).first)
-               {
-                 for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
-                   cell_matrix(i,j) += (((fe_values.shape_value(i,q_point) *
-                                          fe_values.shape_value(j,q_point)) *
-                                         (- om2)
-                                         +
-                                         (fe_values.shape_grad(i,q_point) *
-                                          fe_values.shape_grad(j,q_point)) *
-                                         c2) *
-                                        fe_values.JxW(q_point));
-               }
-           }
-       }
-
-      if (cell->at_boundary())
-       for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-         if (  cell->face(face)->at_boundary() &&
-               (cell->face(face)->boundary_indicator() == 0) )
-           {
-             fe_face_values.reinit (cell, face);
-
-             for (unsigned int i=0; i<dofs_per_cell; ++i)
-               for (unsigned int j=0; j<dofs_per_cell; ++j)
-                 if ((fe.system_to_component_index(i).first != 
-                      fe.system_to_component_index(j).first) &&
-                     fe.has_support_on_face(i, face) &&
-                     fe.has_support_on_face(j, face))
-
-                   for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
-                     cell_matrix(i,j) += ((fe.system_to_component_index(i).first) ? 1 : (-1)) * 
-                                         fe_face_values.shape_value(i,q_point) *
-                                         fe_face_values.shape_value(j,q_point) *
-                                         c * omega *
-                                         fe_face_values.JxW(q_point);
-           }
-
-      cell->get_dof_indices (local_dof_indices);
-
-      for (unsigned int i=0; i<dofs_per_cell; ++i)
-       for (unsigned int j=0; j<dofs_per_cell; ++j)
-         system_matrix.add (   local_dof_indices[i],
-                               local_dof_indices[j],
-                               cell_matrix(i,j));
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+      {
+        if (fe.system_to_component_index(i).first == 
+            fe.system_to_component_index(j).first)
+        {
+          for (unsigned int q_point=0; q_point<n_q_points; ++q_point)
+            cell_matrix(i,j) += (((fe_values.shape_value(i,q_point) *
+                                  fe_values.shape_value(j,q_point)) *
+                                 (- om2)
+                                 +
+                                 (fe_values.shape_grad(i,q_point) *
+                                  fe_values.shape_grad(j,q_point)) *
+                                 c2) *
+                                fe_values.JxW(q_point));
+        }
+      }
     }
 
+    if (cell->at_boundary())
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+        if (cell->face(face)->at_boundary() &&
+            (cell->face(face)->boundary_indicator() == 0) )
+        {
+          fe_face_values.reinit (cell, face);
+
+          for (unsigned int i=0; i<dofs_per_cell; ++i)
+            for (unsigned int j=0; j<dofs_per_cell; ++j)
+              if ((fe.system_to_component_index(i).first != 
+                  fe.system_to_component_index(j).first) &&
+                 fe.has_support_on_face(i, face) &&
+                 fe.has_support_on_face(j, face))
+
+                for (unsigned int q_point=0; q_point<n_face_q_points; ++q_point)
+                  cell_matrix(i,j) += ((fe.system_to_component_index(i).first) ? 1 : (-1)) * 
+                                     fe_face_values.shape_value(i,q_point) *
+                                     fe_face_values.shape_value(j,q_point) *
+                                     c *
+                                     omega *
+                                     fe_face_values.JxW(q_point);
+        }
+
+        cell->get_dof_indices (local_dof_indices);
+
+    for (unsigned int i=0; i<dofs_per_cell; ++i)
+      for (unsigned int j=0; j<dofs_per_cell; ++j)
+        system_matrix.add (local_dof_indices[i],
+                          local_dof_indices[j],
+                          cell_matrix(i,j));
+  }
+
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            1,
                                            DirichletBoundaryValues<dim>(),
                                            boundary_values);
+
   MatrixTools::apply_boundary_values (boundary_values,
                                      system_matrix,
                                      solution,
@@ -338,9 +468,9 @@ void UltrasoundProblem<dim>::assemble_system ()
 
 
 template <int dim>
-void UltrasoundProblem<dim>::solve () 
+void UltrasoundProblem<dim>::solve ()
 {
-  SparseDirectUMFPACK                  A_direct;
+  SparseDirectUMFPACK  A_direct;
 
   A_direct.initialize(system_matrix);
   A_direct.vmult(solution,system_rhs);
@@ -348,56 +478,31 @@ void UltrasoundProblem<dim>::solve ()
 
 
 template <int dim>
-void UltrasoundProblem<dim>::postprocess () 
+void UltrasoundProblem<dim>::output_results () const
 {
-  QMidpoint<dim>       midpoint_rule;
-  FEValues<dim> fe_values(fe, midpoint_rule, update_values);
-
-  std::vector<Vector<double> >         values;
-  values.resize(1);
-  values[0].reinit(2);
+  Postprocessor<dim> pproc;
+  DataOut<dim> data_out;
 
-  typename DoFHandler<dim>::active_cell_iterator
-    cell = dof_handler.begin_active(),
-    endc = dof_handler.end();
-  for (unsigned int    i=0; cell!=endc; ++cell, ++i)
-    {
-      fe_values.reinit(cell);
-               
-      fe_values.get_function_values (solution, values);
+  data_out.attach_dof_handler (dof_handler);
 
-      intensity(i) = sqrt(     values[0](0)*values[0](0) 
-                               + values[0](1)*values[0](1));
-    }
-}
+  prm.enter_subsection("Output parameters");
 
+  const std::string output_file    = prm.get("Output file"),
+                   output_format  = prm.get("Output format");
 
-template <int dim>
-void UltrasoundProblem<dim>::output_results () const
-{
-  const std::string    output_file(prm.get("Output file")), 
-    output_format(prm.get("Output format"));
-       
-  DataOutBase::OutputFormat            format = DataOutBase::parse_output_format(output_format);
+  data_out.parse_parameters(prm);
 
-  const std::string            filename =      output_file + 
-                                               DataOutBase::default_suffix(format);
+  prm.leave_subsection ();
 
-  std::ofstream output (filename.c_str());
+  DataOutBase::OutputFormat  format = DataOutBase::parse_output_format(output_format);
 
-  DataOut<dim> data_out;
-  data_out.parse_parameters(prm);
-  data_out.attach_dof_handler (dof_handler);
+  const std::string filename = output_file +
+                              DataOutBase::default_suffix(format);
 
-  std::vector<std::string> solution_names;
-  solution_names.push_back ("Re_u");
-  solution_names.push_back ("Im_u");
-  data_out.add_data_vector (solution, solution_names);
-       
-  data_out.add_data_vector (intensity, "intensity");
+  std::ofstream output (filename.c_str());
 
+  data_out.add_data_vector (solution, pproc);
   data_out.build_patches ();
-
   data_out.write (output, format);
 }
 
@@ -409,7 +514,6 @@ void UltrasoundProblem<dim>::run ()
   setup_system ();
   assemble_system ();
   solve ();
-  postprocess ();
   output_results ();
 }
 
@@ -417,39 +521,40 @@ void UltrasoundProblem<dim>::run ()
 int main () 
 {
   try
-    {
-      ParameterHandler         prm;
-      ParameterReader          param(prm);
-      param.read_parameters();
-
-      Assert (DIM > 1, ExcNotImplemented());
-               
-      UltrasoundProblem<DIM>   ultrasound_problem (prm);
-      ultrasound_problem.run ();
-    }
+  {
+    ParameterHandler  prm;
+
+    ParameterReader   param(prm);
+    param.read_parameters();
+
+    Assert (DIM > 1, ExcNotImplemented());
+
+    UltrasoundProblem<DIM>  ultrasound_problem (prm);
+    ultrasound_problem.run ();
+  }
   catch (std::exception &exc)
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Exception on processing: " << std::endl
-               << exc.what() << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
+  {
+    std::cerr << std::endl << std::endl
+             << "----------------------------------------------------"
+             << std::endl;
+    std::cerr << "Exception on processing: " << std::endl
+             << exc.what() << std::endl
+             << "Aborting!" << std::endl
+             << "----------------------------------------------------"
+             << std::endl;
+    return 1;
+  }
   catch (...) 
-    {
-      std::cerr << std::endl << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      std::cerr << "Unknown exception!" << std::endl
-               << "Aborting!" << std::endl
-               << "----------------------------------------------------"
-               << std::endl;
-      return 1;
-    }
+  {
+    std::cerr << std::endl << std::endl
+             << "----------------------------------------------------"
+             << std::endl;
+    std::cerr << "Unknown exception!" << std::endl
+             << "Aborting!" << std::endl
+             << "----------------------------------------------------"
+             << std::endl;
+    return 1;
+  }
 
   return 0;
 }
index 9d50174c58a3931bebe53cf9bbce1ddcaaef032f..9f48aef3caf1bc92f48bc08b1b82d541c375cf52 100644 (file)
 # Listing of Parameters
 # ---------------------
-# Distance of the focal point of the lens to the x-axis (or xy-plane in 3D)
-set Focal distance        = 0.3
 
-# Number of global mesh refinement steps applied to initial coarse grid
-set Number of refinements = 5
+subsection Mesh & geometry parameters
+  # Distance of the focal point of the lens to the x-axis (or xy-plane in 3D)
+  set Focal distance        = 0.3
 
-# Name of the output file (without extension)
-set Output file           = solution
-
-# A name for the output format to be used
-set Output format         = vtk
-
-# Number of subdivisions of each mesh cell
-set Subdivisions          = 1
-
-# Wave speed
-set c                     = 1.5e5
-
-# Frequency
-set omega                 = 1.5e7
-
-
-subsection DX output parameters
-  # Output format of vertex coordinates, which is either a text representation
-  # (ascii) or binary floating point values of 32 or 64 bits length
-  set Coordinates format = ascii
-
-  # Output format of data values, which is either a text representation
-  # (ascii) or binary floating point values of 32 or 64 bits length
-  set Data format        = ascii
-
-  # Output format of integer numbers, which is either a text representation
-  # (ascii) or binary integer values of 32 or 64 bits length
-  set Integer format     = ascii
-
-  # A boolean field indicating whether neighborship information between cells
-  # is to be written to the OpenDX output file
-  set Write neighbors    = true
+  # Number of global mesh refinement steps applied to initial coarse grid
+  set Number of refinements = 6
 end
 
 
-subsection Eps output parameters
-  # Angle of the viewing position against the vertical axis
-  set Azimut angle                        = 60
-
-  # Name of a color function used to colorize mesh lines and/or cell
-  # interiors
-  set Color function                      = default
-
-  # Whether the interior of cells shall be shaded
-  set Color shading of interior of cells  = true
-
-  # Whether the mesh lines, or only the surface should be drawn
-  set Draw mesh lines                     = true
-
-  # Whether only the mesh lines, or also the interior of cells should be
-  # plotted. If this flag is false, then one can see through the mesh
-  set Fill interior of cells              = true
+subsection Physical constants
+  # Wave speed
+  set c     = 1.5e5
 
-  # Number of the input vector that is to be used to generate color
-  # information
-  set Index of vector for color           = 0
-
-  # Number of the input vector that is to be used to generate height
-  # information
-  set Index of vector for height          = 0
-
-  # The width in which the postscript renderer is to plot lines
-  set Line widths in eps units            = 0.5
-
-  # Whether width or height should be scaled to match the given size
-  set Scale to width or height            = width
-
-  # Scaling for the z-direction relative to the scaling used in x- and
-  # y-directions
-  set Scaling for z-axis                  = 1
-
-  # The size (width or height) to which the eps output file is to be scaled
-  set Size (width or height) in eps units = 300
-
-  # Angle of the viewing direction against the y-axis
-  set Turn angle                          = 30
+  # Frequency
+  set omega = 5.0e7
 end
 
 
-subsection Povray output parameters
-  # Whether camera and lightling information should be put into an external
-  # file "data.inc" or into the POVRAY input file
-  set Include external file = true
+subsection Output parameters
+  # Name of the output file (without extension)
+  set Output file   = solution
 
-  # Whether POVRAY should use bicubic patches
-  set Use bicubic patches   = false
-
-  # A flag indicating whether POVRAY should use smoothed triangles instead of
-  # the usual ones
-  set Use smooth triangles  = false
-end
-
-
-subsection UCD output parameters
-  # A flag indicating whether a comment should be written to the beginning of
-  # the output file indicating date and time of creation as well as the
-  # creating program
-  set Write preamble = true
+  # A name for the output format to be used
+  set Output format = gmv
 end
-
-

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.