step-29.h

Go to the documentation of this file.

,
 *                               Vector<double> &values) const override
 *     {
 *       Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
 * 
 *       values(0) = 1;
 *       values(1) = 0;
 *     }
 * 
 *     virtual void
 *     vector_value_list(const std::vector<Point<dim>> &points,
 *                       std::vector<Vector<double>> &  value_list) const override
 *     {
 *       Assert(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]);
 *     }
 *   };
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodeParameterReadercodeclass"></a> 
 * <h3>The <code>ParameterReader</code> class</h3>
 * 

 * 
 * The next class is responsible for preparing the ParameterHandler object
 * and reading parameters from an input file.  It includes a function
 * <code>declare_parameters</code> that declares all the necessary
 * parameters and a <code>read_parameters</code> function that is called
 * from outside to initiate the parameter reading process.
 * 
 * @code
 *   class ParameterReader : public Subscriptor
 *   {
 *   public:
 *     ParameterReader(ParameterHandler &);
 *     void read_parameters(const std::string &);
 * 
 *   private:
 *     void              declare_parameters();
 *     ParameterHandler &prm;
 *   };
 * 
 * @endcode
 * 
 * The constructor stores a reference to the ParameterHandler object that is
 * passed to it:
 * 
 * @code
 *   ParameterReader::ParameterReader(ParameterHandler &paramhandler)
 *     : prm(paramhandler)
 *   {}
 * 
 * @endcode
 * 
 * 
 * <a name="codeParameterReaderdeclare_parameterscode"></a> 
 * <h4><code>ParameterReader::declare_parameters</code></h4>
 * 

 * 
 * The <code>declare_parameters</code> function declares all the parameters
 * that our ParameterHandler object will be able to read from input files,
 * along with their types, range conditions and the subsections they appear
 * in. We will wrap all the entries that go into a section in a pair of
 * braces to force the editor to indent them by one level, making it simpler
 * to read which entries together form a section:
 * 
 * @code
 *   void ParameterReader::declare_parameters()
 *   {
 * @endcode
 * 
 * Parameters for mesh and geometry include the number of global
 * refinement steps that are applied to the initial coarse mesh and the
 * focal distance @f$d@f$ of the transducer lens. For the number of refinement
 * steps, we allow integer values in the range @f$[0,\infty)@f$, where the
 * omitted second argument to the Patterns::Integer object denotes the
 * half-open interval.  For the focal distance any number greater than
 * zero is accepted:
 * 
 * @code
 *     prm.enter_subsection("Mesh & geometry parameters");
 *     {
 *       prm.declare_entry("Number of refinements",
 *                         "6",
 *                         Patterns::Integer(0),
 *                         "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");
 *     }
 *     prm.leave_subsection();
 * 
 * @endcode
 * 
 * The next subsection is devoted to the physical parameters appearing in
 * the equation, which are the frequency @f$\omega@f$ and wave speed
 * @f$c@f$. Again, both need to lie in the half-open interval @f$[0,\infty)@f$
 * represented by calling the Patterns::Double class with only the left
 * end-point as argument:
 * 
 * @code
 *     prm.enter_subsection("Physical constants");
 *     {
 *       prm.declare_entry("c", "1.5e5", Patterns::Double(0), "Wave speed");
 * 
 *       prm.declare_entry("omega", "5.0e7", Patterns::Double(0), "Frequency");
 *     }
 *     prm.leave_subsection();
 * 
 * 
 * @endcode
 * 
 * Last but not least we would like to be able to change some properties
 * of the output, like filename and format, through entries in the
 * configuration file, which is the purpose of the last subsection:
 * 
 * @code
 *     prm.enter_subsection("Output parameters");
 *     {
 *       prm.declare_entry("Output filename",
 *                         "solution",
 *                         Patterns::Anything(),
 *                         "Name of the output file (without extension)");
 * 
 * @endcode
 * 
 * Since different output formats may require different parameters for
 * generating output (like for example, postscript output needs
 * viewpoint angles, line widths, colors etc), it would be cumbersome if
 * we had to declare all these parameters by hand for every possible
 * output format supported in the library. Instead, each output format
 * has a <code>FormatFlags::declare_parameters</code> function, which
 * declares all the parameters specific to that format in an own
 * subsection. The following call of
 * DataOutInterface<1>::declare_parameters executes
 * <code>declare_parameters</code> for all available output formats, so
 * that for each format an own subsection will be created with
 * parameters declared for that particular output format. (The actual
 * value of the template parameter in the call, <code>@<1@></code>
 * above, does not matter here: the function does the same work
 * independent of the dimension, but happens to be in a
 * template-parameter-dependent class.)  To find out what parameters
 * there are for which output format, you can either consult the
 * documentation of the DataOutBase class, or simply run this program
 * without a parameter file present. It will then create a file with all
 * declared parameters set to their default values, which can
 * conveniently serve as a starting point for setting the parameters to
 * the values you desire.
 * 
 * @code
 *       DataOutInterface<1>::declare_parameters(prm);
 *     }
 *     prm.leave_subsection();
 *   }
 * 
 * @endcode
 * 
 * 
 * <a name="codeParameterReaderread_parameterscode"></a> 
 * <h4><code>ParameterReader::read_parameters</code></h4>
 * 

 * 
 * This is the main function in the ParameterReader class.  It gets called
 * from outside, first declares all the parameters, and then reads them from
 * the input file whose filename is provided by the caller. After the call
 * to this function is complete, the <code>prm</code> object can be used to
 * retrieve the values of the parameters read in from the file :
 * 
 * @code
 *   void ParameterReader::read_parameters(const std::string &parameter_file)
 *   {
 *     declare_parameters();
 * 
 *     prm.parse_input(parameter_file);
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodeComputeIntensitycodeclass"></a> 
 * <h3>The <code>ComputeIntensity</code> class</h3>
 * 

 * 
 * As mentioned in the introduction, the quantity that we are really after
 * is the spatial distribution of the intensity of the ultrasound wave,
 * which corresponds to @f$|u|=\sqrt{v^2+w^2}@f$. Now we could just be content
 * with having @f$v@f$ and @f$w@f$ in our output, and use a suitable visualization
 * or postprocessing tool to derive @f$|u|@f$ from the solution we
 * computed. However, there is also a way to output data derived from the
 * solution in deal.II, and we are going to make use of this mechanism here.
 * 

 * 
 * So far we have always used the DataOut::add_data_vector function to add
 * vectors containing output data to a DataOut object.  There is a special
 * version of this function that in addition to the data vector has an
 * additional argument of type DataPostprocessor. What happens when this
 * function is used for output is that at each point where output data is to
 * be generated, the DataPostprocessor::evaluate_scalar_field() or
 * DataPostprocessor::evaluate_vector_field() function of the
 * specified DataPostprocessor object is invoked to compute the output
 * quantities from the values, the gradients and the second derivatives of
 * the finite element function represented by the data vector (in the case
 * of face related data, normal vectors are available as well). Hence, this
 * allows us to output any quantity that can locally be derived from the
 * values of the solution and its derivatives.  Of course, the ultrasound
 * intensity @f$|u|@f$ is such a quantity and its computation doesn't even
 * involve any derivatives of @f$v@f$ or @f$w@f$.
 * 

 * 
 * In practice, the DataPostprocessor class only provides an interface to
 * this functionality, and we need to derive our own class from it in order
 * to implement the functions specified by the interface. In the most
 * general case one has to implement several member functions but if the
 * output quantity is a single scalar then some of this boilerplate code can
 * be handled by a more specialized class, DataPostprocessorScalar and we
 * can derive from that one instead. This is what the
 * <code>ComputeIntensity</code> class does:
 * 
 * @code
 *   template <int dim>
 *   class ComputeIntensity : public DataPostprocessorScalar<dim>
 *   {
 *   public:
 *     ComputeIntensity();
 * 
 *     virtual void evaluate_vector_field(
 *       const DataPostprocessorInputs::Vector<dim> &inputs,
 *       std::vector<Vector<double>> &computed_quantities) const override;
 *   };
 * 
 * @endcode
 * 
 * In the constructor, we need to call the constructor of the base class
 * with two arguments. The first denotes the name by which the single scalar
 * quantity computed by this class should be represented in output files. In
 * our case, the postprocessor has @f$|u|@f$ as output, so we use "Intensity".
 *   

 * 
 * The second argument is a set of flags that indicate which data is needed
 * by the postprocessor in order to compute the output quantities.  This can
 * be any subset of update_values, update_gradients and update_hessians
 * (and, in the case of face data, also update_normal_vectors), which are
 * documented in UpdateFlags.  Of course, computation of the derivatives
 * requires additional resources, so only the flags for data that are really
 * needed should be given here, just as we do when we use FEValues objects.
 * In our case, only the function values of @f$v@f$ and @f$w@f$ are needed to
 * compute @f$|u|@f$, so we're good with the update_values flag.
 * 
 * @code
 *   template <int dim>
 *   ComputeIntensity<dim>::ComputeIntensity()
 *     : DataPostprocessorScalar<dim>("Intensity", update_values)
 *   {}
 * 
 * 
 * @endcode
 * 
 * The actual postprocessing happens in the following function. Its input is
 * an object that stores values of the function (which is here vector-valued)
 * representing the data vector given to DataOut::add_data_vector, evaluated
 * at all evaluation points where we generate output, and some tensor objects
 * representing derivatives (that we don't use here since @f$|u|@f$ is computed
 * from just @f$v@f$ and @f$w@f$). The derived quantities are returned in the
 * <code>computed_quantities</code> vector. Remember that this function may
 * only use data for which the respective update flag is specified by
 * <code>get_needed_update_flags</code>. For example, we may not use the
 * derivatives here, since our implementation of
 * <code>get_needed_update_flags</code> requests that only function values
 * are provided.
 * 
 * @code
 *   template <int dim>
 *   void ComputeIntensity<dim>::evaluate_vector_field(
 *     const DataPostprocessorInputs::Vector<dim> &inputs,
 *     std::vector<Vector<double>> &               computed_quantities) const
 *   {
 *     Assert(computed_quantities.size() == inputs.solution_values.size(),
 *            ExcDimensionMismatch(computed_quantities.size(),
 *                                 inputs.solution_values.size()));
 * 
 * @endcode
 * 
 * The computation itself is straightforward: We iterate over each
 * entry in the output vector and compute @f$|u|@f$ from the
 * corresponding values of @f$v@f$ and @f$w@f$. We do this by creating a
 * complex number @f$u@f$ and then calling `std::abs()` on the
 * result. (One may be tempted to call `std::norm()`, but in a
 * historical quirk, the C++ committee decided that `std::norm()`
 * should return the <i>square</i> of the absolute value --
 * thereby not satisfying the properties mathematicians require of
 * something called a "norm".)
 * 
 * @code
 *     for (unsigned int i = 0; i < computed_quantities.size(); i++)
 *       {
 *         Assert(computed_quantities[i].size() == 1,
 *                ExcDimensionMismatch(computed_quantities[i].size(), 1));
 *         Assert(inputs.solution_values[i].size() == 2,
 *                ExcDimensionMismatch(inputs.solution_values[i].size(), 2));
 * 
 *         const std::complex<double> u(inputs.solution_values[i](0),
 *                                      inputs.solution_values[i](1));
 * 
 *         computed_quantities[i](0) = std::abs(u);
 *       }
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="ThecodeUltrasoundProblemcodeclass"></a> 
 * <h3>The <code>UltrasoundProblem</code> class</h3>
 * 

 * 
 * Finally here is the main class of this program.  It's member functions
 * are very similar to the previous examples, in particular @ref step_4 "step-4", and the
 * list of member variables does not contain any major surprises either.
 * The ParameterHandler object that is passed to the constructor is stored
 * as a reference to allow easy access to the parameters from all functions
 * of the class.  Since we are working with vector valued finite elements,
 * the FE object we are using is of type FESystem.
 * 
 * @code
 *   template <int dim>
 *   class UltrasoundProblem
 *   {
 *   public:
 *     UltrasoundProblem(ParameterHandler &);
 *     void run();
 * 
 *   private:
 *     void make_grid();
 *     void setup_system();
 *     void assemble_system();
 *     void solve();
 *     void output_results() const;
 * 
 *     ParameterHandler &prm;
 * 
 *     Triangulation<dim> triangulation;
 *     DoFHandler<dim>    dof_handler;
 *     FESystem<dim>      fe;
 * 
 *     SparsityPattern      sparsity_pattern;
 *     SparseMatrix<double> system_matrix;
 *     Vector<double>       solution, system_rhs;
 *   };
 * 
 * 
 * 
 * @endcode
 * 
 * The constructor takes the ParameterHandler object and stores it in a
 * reference. It also initializes the DoF-Handler and the finite element
 * system, which consists of two copies of the scalar Q1 field, one for @f$v@f$
 * and one for @f$w@f$:
 * 
 * @code
 *   template <int dim>
 *   UltrasoundProblem<dim>::UltrasoundProblem(ParameterHandler &param)
 *     : prm(param)
 *     , dof_handler(triangulation)
 *     , fe(FE_Q<dim>(1), 2)
 *   {}
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemmake_gridcode"></a> 
 * <h4><code>UltrasoundProblem::make_grid</code></h4>
 * 

 * 
 * Here we setup the grid for our domain.  As mentioned in the exposition,
 * the geometry is just a unit square (in 2d) with the part of the boundary
 * that represents the transducer lens replaced by a sector of a circle.
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::make_grid()
 *   {
 * @endcode
 * 
 * First we generate some logging output and start a timer so we can
 * compute execution time when this function is done:
 * 
 * @code
 *     deallog << "Generating grid... ";
 *     Timer timer;
 * 
 * @endcode
 * 
 * Then we query the values for the focal distance of the transducer lens
 * and the number of mesh refinement steps from our ParameterHandler
 * object:
 * 
 * @code
 *     prm.enter_subsection("Mesh & geometry parameters");
 * 
 *     const double       focal_distance = prm.get_double("Focal distance");
 *     const unsigned int n_refinements = prm.get_integer("Number of refinements");
 * 
 *     prm.leave_subsection();
 * 
 * @endcode
 * 
 * Next, two points are defined for position and focal point of the
 * transducer lens, which is the center of the circle whose segment will
 * form the transducer part of the boundary. Notice that this is the only
 * point in the program where things are slightly different in 2D and 3D.
 * Even though this tutorial only deals with the 2D case, the necessary
 * additions to make this program functional in 3D are so minimal that we
 * opt for including them:
 * 
 * @code
 *     const Point<dim> transducer =
 *       (dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
 *     const Point<dim> focal_point = (dim == 2) ?
 *                                      Point<dim>(0.5, focal_distance) :
 *                                      Point<dim>(0.5, 0.5, focal_distance);
 * 
 * 
 * @endcode
 * 
 * As initial coarse grid we take a simple unit square with 5 subdivisions
 * in each direction. The number of subdivisions is chosen so that the
 * line segment @f$[0.4,0.6]@f$ that we want to designate as the transducer
 * boundary is spanned by a single face. Then we step through all cells to
 * find the faces where the transducer is to be located, which in fact is
 * just the single edge from 0.4 to 0.6 on the x-axis. This is where we
 * want the refinements to be made according to a circle shaped boundary,
 * so we mark this edge with a different manifold indicator. Since we will
 * Dirichlet boundary conditions on the transducer, we also change its
 * boundary indicator.
 * 
 * @code
 *     GridGenerator::subdivided_hyper_cube(triangulation, 5, 0, 1);
 * 
 *     for (auto &cell : triangulation.cell_iterators())
 *       for (const auto &face : cell->face_iterators())
 *         if (face->at_boundary() &&
 *             ((face->center() - transducer).norm_square() < 0.01))
 *           {
 *             face->set_boundary_id(1);
 *             face->set_manifold_id(1);
 *           }
 * @endcode
 * 
 * For the circle part of the transducer lens, a SphericalManifold object
 * is used (which, of course, in 2D just represents a circle), with center
 * computed as above.
 * 
 * @code
 *     triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
 * 
 * @endcode
 * 
 * Now global refinement is executed. Cells near the transducer location
 * will be automatically refined according to the circle shaped boundary
 * of the transducer lens:
 * 
 * @code
 *     triangulation.refine_global(n_refinements);
 * 
 * @endcode
 * 
 * Lastly, we generate some more logging output. We stop the timer and
 * query the number of CPU seconds elapsed since the beginning of the
 * function:
 * 
 * @code
 *     timer.stop();
 *     deallog << "done (" << timer.cpu_time() << "s)" << std::endl;
 * 
 *     deallog << "  Number of active cells:  " << triangulation.n_active_cells()
 *             << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemsetup_systemcode"></a> 
 * <h4><code>UltrasoundProblem::setup_system</code></h4>
 *   

 * 
 * Initialization of the system matrix, sparsity patterns and vectors are
 * the same as in previous examples and therefore do not need further
 * comment. As in the previous function, we also output the run time of what
 * we do here:
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::setup_system()
 *   {
 *     deallog << "Setting up system... ";
 *     Timer timer;
 * 
 *     dof_handler.distribute_dofs(fe);
 * 
 *     DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
 *     DoFTools::make_sparsity_pattern(dof_handler, dsp);
 *     sparsity_pattern.copy_from(dsp);
 * 
 *     system_matrix.reinit(sparsity_pattern);
 *     system_rhs.reinit(dof_handler.n_dofs());
 *     solution.reinit(dof_handler.n_dofs());
 * 
 *     timer.stop();
 *     deallog << "done (" << timer.cpu_time() << "s)" << std::endl;
 * 
 *     deallog << "  Number of degrees of freedom: " << dof_handler.n_dofs()
 *             << std::endl;
 *   }
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemassemble_systemcode"></a> 
 * <h4><code>UltrasoundProblem::assemble_system</code></h4>
 * 

 * 
 * As before, this function takes care of assembling the system matrix and
 * right hand side vector:
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::assemble_system()
 *   {
 *     deallog << "Assembling system matrix... ";
 *     Timer timer;
 * 
 * @endcode
 * 
 * First we query wavespeed and frequency from the ParameterHandler object
 * and store them in local variables, as they will be used frequently
 * throughout this function.
 * 

 * 
 * 
 * @code
 *     prm.enter_subsection("Physical constants");
 * 
 *     const double omega = prm.get_double("omega"), c = prm.get_double("c");
 * 
 *     prm.leave_subsection();
 * 
 * @endcode
 * 
 * As usual, for computing integrals ordinary Gauss quadrature rule is
 * used. Since our bilinear form involves boundary integrals on
 * @f$\Gamma_2@f$, we also need a quadrature rule for surface integration on
 * the faces, which are @f$dim-1@f$ dimensional:
 * 
 * @code
 *     QGauss<dim>     quadrature_formula(fe.degree + 1);
 *     QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
 * 
 *     const unsigned int n_q_points      = quadrature_formula.size(),
 *                        n_face_q_points = face_quadrature_formula.size(),
 *                        dofs_per_cell   = fe.dofs_per_cell;
 * 
 * @endcode
 * 
 * The FEValues objects will evaluate the shape functions for us.  For the
 * part of the bilinear form that involves integration on @f$\Omega@f$, we'll
 * need the values and gradients of the shape functions, and of course the
 * quadrature weights.  For the terms involving the boundary integrals,
 * only shape function values and the quadrature weights are necessary.
 * 
 * @code
 *     FEValues<dim> fe_values(fe,
 *                             quadrature_formula,
 *                             update_values | update_gradients |
 *                               update_JxW_values);
 * 
 *     FEFaceValues<dim> fe_face_values(fe,
 *                                      face_quadrature_formula,
 *                                      update_values | update_JxW_values);
 * 
 * @endcode
 * 
 * As usual, the system matrix is assembled cell by cell, and we need a
 * matrix for storing the local cell contributions as well as an index
 * vector to transfer the cell contributions to the appropriate location
 * in the global system matrix after.
 * 
 * @code
 *     FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *     std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 * 
 *     for (const auto &cell : dof_handler.active_cell_iterators())
 *       {
 * @endcode
 * 
 * On each cell, we first need to reset the local contribution matrix
 * and request the FEValues object to compute the shape functions for
 * the current cell:
 * 
 * @code
 *         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)
 *               {
 * @endcode
 * 
 * At this point, it is important to keep in mind that we are
 * dealing with a finite element system with two
 * components. Due to the way we constructed this FESystem,
 * namely as the Cartesian product of two scalar finite
 * element fields, each shape function has only a single
 * nonzero component (they are, in deal.II lingo, @ref
 * GlossPrimitive "primitive").  Hence, each shape function
 * can be viewed as one of the @f$\phi@f$'s or @f$\psi@f$'s from the
 * introduction, and similarly the corresponding degrees of
 * freedom can be attributed to either @f$\alpha@f$ or @f$\beta@f$.
 * As we iterate through all the degrees of freedom on the
 * current cell however, they do not come in any particular
 * order, and so we cannot decide right away whether the DoFs
 * with index @f$i@f$ and @f$j@f$ belong to the real or imaginary part
 * of our solution.  On the other hand, if you look at the
 * form of the system matrix in the introduction, this
 * distinction is crucial since it will determine to which
 * block in the system matrix the contribution of the current
 * pair of DoFs will go and hence which quantity we need to
 * compute from the given two shape functions.  Fortunately,
 * the FESystem object can provide us with this information,
 * namely it has a function
 * FESystem::system_to_component_index, that for each local
 * DoF index returns a pair of integers of which the first
 * indicates to which component of the system the DoF
 * belongs. The second integer of the pair indicates which
 * index the DoF has in the scalar base finite element field,
 * but this information is not relevant here. If you want to
 * know more about this function and the underlying scheme
 * behind primitive vector valued elements, take a look at
 * @ref step_8 "step-8" or the @ref vector_valued module, where these topics
 * are explained in depth.
 * 
 * @code
 *                 if (fe.system_to_component_index(i).first ==
 *                     fe.system_to_component_index(j).first)
 *                   {
 * @endcode
 * 
 * If both DoFs @f$i@f$ and @f$j@f$ belong to same component,
 * i.e. their shape functions are both @f$\phi@f$'s or both
 * @f$\psi@f$'s, the contribution will end up in one of the
 * diagonal blocks in our system matrix, and since the
 * corresponding entries are computed by the same formula,
 * we do not bother if they actually are @f$\phi@f$ or @f$\psi@f$
 * shape functions. We can simply compute the entry by
 * iterating over all quadrature points and adding up
 * their contributions, where values and gradients of the
 * shape functions are supplied by our FEValues object.
 * 

 * 
 * 
 * @code
 *                     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)) *
 *                             (-omega * omega) +
 *                           (fe_values.shape_grad(i, q_point) *
 *                            fe_values.shape_grad(j, q_point)) *
 *                             c * c) *
 *                          fe_values.JxW(q_point));
 * 
 * @endcode
 * 
 * You might think that we would have to specify which
 * component of the shape function we'd like to evaluate
 * when requesting shape function values or gradients from
 * the FEValues object. However, as the shape functions
 * are primitive, they have only one nonzero component,
 * and the FEValues class is smart enough to figure out
 * that we are definitely interested in this one nonzero
 * component.
 * 
 * @code
 *                   }
 *               }
 *           }
 * 
 * 
 * @endcode
 * 
 * We also have to add contributions due to boundary terms. To this
 * end, we loop over all faces of the current cell and see if first it
 * is at the boundary, and second has the correct boundary indicator
 * associated with @f$\Gamma_2@f$, the part of the boundary where we have
 * absorbing boundary conditions:
 * 
 * @code
 *         for (unsigned int face_no : GeometryInfo<dim>::face_indices())
 *           if (cell->face(face_no)->at_boundary() &&
 *               (cell->face(face_no)->boundary_id() == 0))
 *             {
 * @endcode
 * 
 * These faces will certainly contribute to the off-diagonal
 * blocks of the system matrix, so we ask the FEFaceValues
 * object to provide us with the shape function values on this
 * face:
 * 
 * @code
 *               fe_face_values.reinit(cell, face_no);
 * 
 * 
 * @endcode
 * 
 * Next, we loop through all DoFs of the current cell to find
 * pairs that belong to different components and both have
 * support on the current face_no:
 * 
 * @code
 *               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_no) &&
 *                       fe.has_support_on_face(j, face_no))
 * @endcode
 * 
 * The check whether shape functions have support on a
 * face is not strictly necessary: if we don't check for
 * it we would simply add up terms to the local cell
 * matrix that happen to be zero because at least one of
 * the shape functions happens to be zero. However, we can
 * save that work by adding the checks above.
 * 

 * 
 * In either case, these DoFs will contribute to the
 * boundary integrals in the off-diagonal blocks of the
 * system matrix. To compute the integral, we loop over
 * all the quadrature points on the face and sum up the
 * contribution weighted with the quadrature weights that
 * the face quadrature rule provides.  In contrast to the
 * entries on the diagonal blocks, here it does matter
 * which one of the shape functions is a @f$\psi@f$ and which
 * one is a @f$\phi@f$, since that will determine the sign of
 * the entry.  We account for this by a simple conditional
 * statement that determines the correct sign. Since we
 * already checked that DoF @f$i@f$ and @f$j@f$ belong to
 * different components, it suffices here to test for one
 * of them to which component it belongs.
 * 
 * @code
 *                     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 == 0) ? -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);
 *             }
 * 
 * @endcode
 * 
 * Now we are done with this cell and have to transfer its
 * contributions from the local to the global system matrix. To this
 * end, we first get a list of the global indices of the this cells
 * DoFs...
 * 
 * @code
 *         cell->get_dof_indices(local_dof_indices);
 * 
 * 
 * @endcode
 * 
 * ...and then add the entries to the system matrix one by one:
 * 
 * @code
 *         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));
 *       }
 * 
 * 
 * @endcode
 * 
 * The only thing left are the Dirichlet boundary values on @f$\Gamma_1@f$,
 * which is characterized by the boundary indicator 1. The Dirichlet
 * values are provided by the <code>DirichletBoundaryValues</code> class
 * we defined above:
 * 
 * @code
 *     std::map<types::global_dof_index, double> boundary_values;
 *     VectorTools::interpolate_boundary_values(dof_handler,
 *                                              1,
 *                                              DirichletBoundaryValues<dim>(),
 *                                              boundary_values);
 * 
 *     MatrixTools::apply_boundary_values(boundary_values,
 *                                        system_matrix,
 *                                        solution,
 *                                        system_rhs);
 * 
 *     timer.stop();
 *     deallog << "done (" << timer.cpu_time() << "s)" << std::endl;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemsolvecode"></a> 
 * <h4><code>UltrasoundProblem::solve</code></h4>
 * 

 * 
 * As already mentioned in the introduction, the system matrix is neither
 * symmetric nor definite, and so it is not quite obvious how to come up
 * with an iterative solver and a preconditioner that do a good job on this
 * matrix.  We chose instead to go a different way and solve the linear
 * system with the sparse LU decomposition provided by UMFPACK. This is
 * often a good first choice for 2D problems and works reasonably well even
 * for a large number of DoFs.  The deal.II interface to UMFPACK is given by
 * the SparseDirectUMFPACK class, which is very easy to use and allows us to
 * solve our linear system with just 3 lines of code.
 * 

 * 
 * Note again that for compiling this example program, you need to have the
 * deal.II library built with UMFPACK support.
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::solve()
 *   {
 *     deallog << "Solving linear system... ";
 *     Timer timer;
 * 
 * @endcode
 * 
 * The code to solve the linear system is short: First, we allocate an
 * object of the right type. The following <code>initialize</code> call
 * provides the matrix that we would like to invert to the
 * SparseDirectUMFPACK object, and at the same time kicks off the
 * LU-decomposition. Hence, this is also the point where most of the
 * computational work in this program happens.
 * 
 * @code
 *     SparseDirectUMFPACK A_direct;
 *     A_direct.initialize(system_matrix);
 * 
 * @endcode
 * 
 * After the decomposition, we can use <code>A_direct</code> like a matrix
 * representing the inverse of our system matrix, so to compute the
 * solution we just have to multiply with the right hand side vector:
 * 
 * @code
 *     A_direct.vmult(solution, system_rhs);
 * 
 *     timer.stop();
 *     deallog << "done (" << timer.cpu_time() << "s)" << std::endl;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemoutput_resultscode"></a> 
 * <h4><code>UltrasoundProblem::output_results</code></h4>
 * 

 * 
 * Here we output our solution @f$v@f$ and @f$w@f$ as well as the derived quantity
 * @f$|u|@f$ in the format specified in the parameter file. Most of the work for
 * deriving @f$|u|@f$ from @f$v@f$ and @f$w@f$ was already done in the implementation of
 * the <code>ComputeIntensity</code> class, so that the output routine is
 * rather straightforward and very similar to what is done in the previous
 * tutorials.
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::output_results() const
 *   {
 *     deallog << "Generating output... ";
 *     Timer timer;
 * 
 * @endcode
 * 
 * Define objects of our <code>ComputeIntensity</code> class and a DataOut
 * object:
 * 
 * @code
 *     ComputeIntensity<dim> intensities;
 *     DataOut<dim>          data_out;
 * 
 *     data_out.attach_dof_handler(dof_handler);
 * 
 * @endcode
 * 
 * Next we query the output-related parameters from the ParameterHandler.
 * The DataOut::parse_parameters call acts as a counterpart to the
 * DataOutInterface<1>::declare_parameters call in
 * <code>ParameterReader::declare_parameters</code>. It collects all the
 * output format related parameters from the ParameterHandler and sets the
 * corresponding properties of the DataOut object accordingly.
 * 
 * @code
 *     prm.enter_subsection("Output parameters");
 * 
 *     const std::string output_filename = prm.get("Output filename");
 *     data_out.parse_parameters(prm);
 * 
 *     prm.leave_subsection();
 * 
 * @endcode
 * 
 * Now we put together the filename from the base name provided by the
 * ParameterHandler and the suffix which is provided by the DataOut class
 * (the default suffix is set to the right type that matches the one set
 * in the .prm file through parse_parameters()):
 * 
 * @code
 *     const std::string filename = output_filename + data_out.default_suffix();
 * 
 *     std::ofstream output(filename);
 * 
 * @endcode
 * 
 * The solution vectors @f$v@f$ and @f$w@f$ are added to the DataOut object in the
 * usual way:
 * 
 * @code
 *     std::vector<std::string> solution_names;
 *     solution_names.emplace_back("Re_u");
 *     solution_names.emplace_back("Im_u");
 * 
 *     data_out.add_data_vector(solution, solution_names);
 * 
 * @endcode
 * 
 * For the intensity, we just call <code>add_data_vector</code> again, but
 * this with our <code>ComputeIntensity</code> object as the second
 * argument, which effectively adds @f$|u|@f$ to the output data:
 * 
 * @code
 *     data_out.add_data_vector(solution, intensities);
 * 
 * @endcode
 * 
 * The last steps are as before. Note that the actual output format is now
 * determined by what is stated in the input file, i.e. one can change the
 * output format without having to re-compile this program:
 * 
 * @code
 *     data_out.build_patches();
 *     data_out.write(output);
 * 
 *     timer.stop();
 *     deallog << "done (" << timer.cpu_time() << "s)" << std::endl;
 *   }
 * 
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="codeUltrasoundProblemruncode"></a> 
 * <h4><code>UltrasoundProblem::run</code></h4>
 * 

 * 
 * Here we simply execute our functions one after the other:
 * 
 * @code
 *   template <int dim>
 *   void UltrasoundProblem<dim>::run()
 *   {
 *     make_grid();
 *     setup_system();
 *     assemble_system();
 *     solve();
 *     output_results();
 *   }
 * } // namespace Step29
 * 
 * 
 * @endcode
 * 
 * 
 * <a name="Thecodemaincodefunction"></a> 
 * <h4>The <code>main</code> function</h4>
 * 

 * 
 * Finally the <code>main</code> function of the program. It has the same
 * structure as in almost all of the other tutorial programs. The only
 * exception is that we define ParameterHandler and
 * <code>ParameterReader</code> objects, and let the latter read in the
 * parameter values from a textfile called <code>@ref step_29 "step-29".prm</code>. The
 * values so read are then handed over to an instance of the UltrasoundProblem
 * class:
 * 
 * @code
 * int main()
 * {
 *   try
 *     {
 *       using namespace dealii;
 *       using namespace Step29;
 * 
 *       deallog.depth_console(5);
 * 
 *       ParameterHandler prm;
 *       ParameterReader  param(prm);
 *       param.read_parameters("step-29.prm");
 * 
 *       UltrasoundProblem<2> 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;
 *     }
 *   catch (...)
 *     {
 *       std::cerr << std::endl
 *                 << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       std::cerr << "Unknown exception!" << std::endl
 *                 << "Aborting!" << std::endl
 *                 << "----------------------------------------------------"
 *                 << std::endl;
 *       return 1;
 *     }
 *   return 0;
 * }
 * @endcode
<a name="Results"></a>
<a name="Results"></a><h1>Results</h1>


The current program reads its run-time parameters from an input file
called <code>step-29.prm</code> that looks like this:
@code
subsection Mesh & geometry parameters
  # Distance of the focal point of the lens to the x-axis
  set Focal distance        = 0.3

  # Number of global mesh refinement steps applied to initial coarse grid
  set Number of refinements = 5
end


subsection Physical constants
  # Wave speed
  set c     = 1.5e5

  # Frequency
  set omega = 3.0e7
end


subsection Output parameters
  # Name of the output file (without extension)
  set Output file   = solution

  # A name for the output format to be used
  set Output format = vtu
end
@endcode

As can be seen, we set
@f$d=0.3@f$, which amounts to a focus of the transducer lens
at @f$x=0.5@f$, @f$y=0.3@f$. The coarse mesh is refined 5 times,
resulting in 160x160 cells, and the output is written in vtu
format. The parameter reader understands many more parameters
pertaining in particular to the generation of output, but we
need none of these parameters here and therefore stick with
their default values.

Here's the console output of the program in debug mode:

@code
> make run
[ 66%] Built target @ref step_29 "step-29"
[100%] Run @ref step_29 "step-29" with Debug configuration
DEAL::Generating grid... done (0.832682s)
DEAL::  Number of active cells:  25600
DEAL::Setting up system... done (0.614761s)
DEAL::  Number of degrees of freedom: 51842
DEAL::Assembling system matrix... done (2.82536s)
DEAL::Solving linear system... done (2.27971s)
DEAL::Generating output... done (1.84418s)
[100%] Built target run
@endcode

(Of course, execution times will differ if you run the program
locally.) The fact that most of the time is spent on assembling
the system matrix and generating output is due to the many assertions
that need to be checked in debug mode. In release mode these parts
of the program run much faster whereas solving the linear system is
hardly sped up at all:

@code
> make run
[ 66%] Built target @ref step_29 "step-29"
Scanning dependencies of target run
[100%] Run @ref step_29 "step-29" with Release configuration
DEAL::Generating grid... done (0.0144960s)
DEAL::  Number of active cells:  25600
DEAL::Setting up system... done (0.0356880s)
DEAL::  Number of degrees of freedom: 51842
DEAL::Assembling system matrix... done (0.0436570s)
DEAL::Solving linear system... done (1.54733s)
DEAL::Generating output... done (0.720528s)
[100%] Built target run
@endcode

The graphical output of the program looks as follows:


<table align="center" class="doxtable">
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-29.v.png" alt="v = Re(u)">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-29.w.png" alt="w = Im(u)">
    </td>
  </tr>
  <tr>
    <td colspan="2">
      <img src="https://www.dealii.org/images/steps/developer/step-29.intensity.png" alt="|u|">
    </td>
  </tr>
</table>

The first two pictures show the real and imaginary parts of
@f$u@f$, whereas the last shows the intensity @f$|u|@f$. One can clearly
see that the intensity is focused around the focal point of the
lens (0.5, 0.3), and that the focus
is rather sharp in @f$x@f$-direction but more blurred in @f$y@f$-direction, which is a
consequence of the geometry of the focusing lens, its finite aperture,
and the wave nature of the problem.

Because colorful graphics are always fun, and to stress the focusing
effects some more, here is another set of images highlighting how well
the intensity is actually focused in @f$x@f$-direction:

<table align="center" class="doxtable">
  <tr>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-29.surface.png" alt="|u|">
    </td>
    <td>
      <img src="https://www.dealii.org/images/steps/developer/step-29.contours.png" alt="|u|">
    </td>
  </tr>
</table>


As a final note, the structure of the program makes it easy to
determine which parts of the program scale nicely as the mesh is
refined and which parts don't. Here are the run times for 5, 6, and 7
global refinements:

@code
> make run
[ 66%] Built target @ref step_29 "step-29"
[100%] Run @ref step_29 "step-29" with Release configuration
DEAL::Generating grid... done (0.0135260s)
DEAL::  Number of active cells:  25600
DEAL::Setting up system... done (0.0213910s)
DEAL::  Number of degrees of freedom: 51842
DEAL::Assembling system matrix... done (0.0414300s)
DEAL::Solving linear system... done (1.56621s)
DEAL::Generating output... done (0.729605s)
[100%] Built target run

> make run
[ 66%] Built target @ref step_29 "step-29"
[100%] Run @ref step_29 "step-29" with Release configuration
DEAL::Generating grid... done (0.0668490s)
DEAL::  Number of active cells:  102400
DEAL::Setting up system... done (0.109694s)
DEAL::  Number of degrees of freedom: 206082
DEAL::Assembling system matrix... done (0.160784s)
DEAL::Solving linear system... done (7.86577s)
DEAL::Generating output... done (2.89320s)
[100%] Built target run

> make run
[ 66%] Built target @ref step_29 "step-29"
[100%] Run @ref step_29 "step-29" with Release configuration
DEAL::Generating grid... done (0.293154s)
DEAL::  Number of active cells:  409600
DEAL::Setting up system... done (0.491301s)
DEAL::  Number of degrees of freedom: 821762
DEAL::Assembling system matrix... done (0.605386s)
DEAL::Solving linear system... done (45.1989s)
DEAL::Generating output... done (11.2292s)
[100%] Built target run
@endcode

Each time we refine the mesh once, so the number of cells and degrees
of freedom roughly quadruples from each step to the next. As can be seen,
generating the grid, setting up degrees of freedom, assembling the
linear system, and generating output scale pretty closely to linear,
whereas solving the linear system is an operation that requires 8
times more time each time the number of degrees of freedom is
increased by a factor of 4, i.e. it is @f${\cal O}(N^{3/2})@f$. This can
be explained by the fact that (using optimal ordering) the
bandwidth of a finite element matrix is @f$B={\cal O}(N^{(dim-1)/dim})@f$,
and the effort to solve a banded linear system using LU decomposition
is @f${\cal O}(BN)@f$. This also explains why the program does run in 3d
as well (after changing the dimension on the
<code>UltrasoundProblem</code> object), but scales very badly and
takes extraordinary patience before it finishes solving the linear
system on a mesh with appreciable resolution, even though all the
other parts of the program scale very nicely.



<a name="extensions"></a>
<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>


An obvious possible extension for this program is to run it in 3d
&mdash; after all, the world around us is three-dimensional, and
ultrasound beams propagate in three-dimensional media. You can try
this by simply changing the template parameter of the principal class
in <code>main()</code> and running it. This won't get you very far,
though: certainly not if you do 5 global refinement steps as set in
the parameter file. You'll simply run out of memory as both the mesh
(with its @f$(2^5)^3 \cdot 5^3=2^{15}\cdot 125 \approx 4\cdot 10^6@f$ cells)
and in particular the sparse direct solver take too much memory. You
can solve with 3 global refinement steps, however, if you have a bit
of time: in early 2011, the direct solve takes about half an
hour. What you'll notice, however, is that the solution is completely
wrong: the mesh size is simply not small enough to resolve the
solution's waves accurately, and you can see this in plots of the
solution. Consequently, this is one of the cases where adaptivity is
indispensable if you don't just want to throw a bigger (presumably
%parallel) machine at the problem.
 *
 *
<a name="PlainProg"></a>
<h1> The plain program</h1>
@include "step-29.cc"
*/