]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Some small changes about step-29
authorMoritz Allmaras <allmaras@math.tamu.edu>
Wed, 31 Oct 2007 16:35:40 +0000 (16:35 +0000)
committerMoritz Allmaras <allmaras@math.tamu.edu>
Wed, 31 Oct 2007 16:35:40 +0000 (16:35 +0000)
git-svn-id: https://svn.dealii.org/trunk@15416 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/doc/authors.html
deal.II/doc/news/changes.h
deal.II/doc/news/news.html
deal.II/examples/step-29/doc/results.dox
deal.II/examples/step-29/step-29.cc

index cf24a1e525513b196a13a81b290f29211078e216..3214bfa8dab1972c31544bbcd6c50723fe23060e 100644 (file)
@@ -53,6 +53,9 @@ alphabetical order):
 </p>
 
 <ul>
+<li><em>Moritz Allmaras:</em>
+    Step-29 tutorial program.
+
 <li><em>Michael Anderson:</em>
     Linear complexity grid reordering algorithm in 2d and 3d.
 
index f51ed6d4638637dc5ba43339d6df3965f3572714..7189c9fdcf948045f953c7ac7fe8c4ead6e2a9cf 100644 (file)
@@ -84,6 +84,12 @@ inconvenience this causes.
 <h3>General</h3>
 
 <ol>
+  <li> <p>New: step-27 demonstrates how problems involving complex numbers 
+  can be implemented by viewing real and imaginary parts of a complex-valued 
+  solution as the two components of a vector-valued function. 
+  <br>
+  (Moritz Allmaras, 2007/10/31)
+  </p>
 
   <li> <p>Improved: On Mac OS X, the operating system provides for
   "frameworks", which are essentially collections of shared libraries.
@@ -95,7 +101,8 @@ inconvenience this causes.
   or vectorized libraries if they are available.
   <br>
   (Eh Tan, WB 2007/10/22)
-  </p></li>
+  </p>
+
 </ol>
 
 
index 85528d42308802be0852032712c3904b6eee1a61..306ae50e1e0cb0ac22124952af76e011f22df7d4 100644 (file)
     <a href="../doxygen/deal.II/changes_after_6.0.html">here.</a></p>
 
     <dl>
+      <dt>
+          <strong>2007/10/31: step-29 tutorial added</strong>
+      </dt>
+      <dd>
+          The step-29 tutorial program explains how problems involving complex 
+          numbers can be treated in <acronym>deal.II</acronym>. In the example, 
+          the complex-valued solution to a Helmholtz equation is split into 
+          real and imaginary parts and a vector-valued finite element field 
+          is used to discretize them seperately. The program also demonstrates 
+          how quantities derived from the finite element solution can be directly 
+          included in the output by using the new <code>DataPostprocessor</code>
+          class. 
+      </dd>
+
       <dt>
           <strong style="color:red">2007/09/07: Version 6.0 released</strong>
       </dt>
index f59dbabf430eb192c534c77dda75ecff7c8f932d..01a5febdb29ff1a0da6b81b92bba4d1c705f88d4 100644 (file)
@@ -148,7 +148,7 @@ DEAL::Generating output... done (15.1489s)
 @endcode
 
 Each time we refine the mesh once, so the number of cells and degrees
-of freedom roughlt quadruples from each step to the next. As can be seen, 
+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
index ed9f27c8a2c77a676efc8349a984bc2be7083796..f9c34d7b8fa88266c86d76633c829e7411568ae5 100644 (file)
 
 
 
-                               // @sect3{Include files}
+                                // @sect3{Include files}
 
-                               // The following header files are unchanged 
-                               // from step-7 and have been discussed before:
+                                // The following header files are unchanged 
+                                // from step-7 and have been discussed before:
 
 #include <base/quadrature_lib.h>
 #include <base/function.h>
 
 #include <fstream>
 
-                               // This header file contains the necessary 
-                               // declarations for the ParameterHandler class
-                               // that we will use to read our parameters from 
-                               // a configuration file: 
+                                // This header file contains the
+                                // necessary declarations for the
+                                // ParameterHandler class that we
+                                // will use to read our parameters
+                                // from a configuration file:
 #include <base/parameter_handler.h>
 
-                               // For solving the linear system, we'll use 
-                               // the sparse LU-decomposition provided by UMFPACK, 
-                               // for which the following header file is needed. 
-                               // Note that in order to compile this tutorial program, 
-                               // the deal.II-library needs to be 
-                               // built with UMFPACK support, which can be most easily 
-                               // achieved by giving the <code> --with-umfpack</code>
-                               // switch when configuring the library:
+                                // For solving the linear system,
+                                // we'll use the sparse
+                                // LU-decomposition provided by
+                                // UMFPACK, for which the following
+                                // header file is needed.  Note that
+                                // in order to compile this tutorial
+                                // program, the deal.II-library needs
+                                // to be built with UMFPACK support,
+                                // which can be most easily achieved
+                                // by giving the <code>
+                                // --with-umfpack</code> switch when
+                                // configuring the library:
 #include <lac/sparse_direct.h>
 
-                               // The FESystem class allows us to stack
-                               // several FE-objects to one compound, vector-valued finite
-                               // element field. The necessary declarations for this class
-                               // are provided in this header file:
+                                // The FESystem class allows us to
+                                // stack several FE-objects to one
+                                // compound, vector-valued finite
+                                // element field. The necessary
+                                // declarations for this class are
+                                // provided in this header file:
 #include <fe/fe_system.h>
 
-                               // Finally, include the header file that declares the
-                               // Timer class that we will use to determine how
-                               // much time each of the operations of our program
-                               // takes:
+                                // Finally, include the header file
+                                // that declares the Timer class that
+                                // we will use to determine how much
+                                // time each of the operations of our
+                                // program takes:
 #include <base/timer.h>
 
-                               // As the last step at the beginning of this program,
-                               // we make everything that is in the
-                               // deal.II namespace globally
-                               // available, without the need to
-                               // prefix everything with
-                               // <code>dealii</code><code>::</code>:
+                                // As the last step at the beginning
+                                // of this program, we make
+                                // everything that is in the deal.II
+                                // namespace globally available,
+                                // without the need to prefix
+                                // everything with
+                                // <code>dealii</code><code>::</code>:
 using namespace dealii;
 
 
-                               // @sect3{The <code>DirichletBoundaryValues</code> class}
-
-                               // First we define a class for the function representing 
-                               // the Dirichlet boundary values. This has been done many times before
-                               // and therefore does not need much explanation.
-                               //
-                               // Since there are two values $v$ and
-                               // $w$ that need to be prescribed at
-                               // the boundary, we have to tell the
-                               // base class that this is a
-                               // vector-valued function with two
-                               // components, and the
-                               // <code>vector_value</code> function
-                               // and its cousin
-                               // <code>vector_value_list</code> must
-                               // return vectors with two entries. In
-                               // our case the function is very
-                               // simple, it just returns 1 for the
-                               // real part $v$ and 0 for the
-                               // imaginary part $w$ regardless of
-                               // the point where it is evaluated.
+                                // @sect3{The <code>DirichletBoundaryValues</code> class}
+
+                                // First we define a class for the
+                                // function representing the
+                                // Dirichlet boundary values. This
+                                // has been done many times before
+                                // and therefore does not need much
+                                // explanation.
+                                //
+                                // Since there are two values $v$ and
+                                // $w$ that need to be prescribed at
+                                // the boundary, we have to tell the
+                                // base class that this is a
+                                // vector-valued function with two
+                                // components, and the
+                                // <code>vector_value</code> function
+                                // and its cousin
+                                // <code>vector_value_list</code> must
+                                // return vectors with two entries. In
+                                // our case the function is very
+                                // simple, it just returns 1 for the
+                                // real part $v$ and 0 for the
+                                // imaginary part $w$ regardless of
+                                // the point where it is evaluated.
 template <int dim>
 class DirichletBoundaryValues : public Function<dim>
 {
@@ -138,16 +150,20 @@ void DirichletBoundaryValues<dim>::vector_value_list (const std::vector<Point<di
     DirichletBoundaryValues<dim>::vector_value (points[p], value_list[p]);
 }
 
-                               // @sect3{The <code>ParameterReader</code> class}
-
-                               // 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. 
+                                // @sect3{The <code>ParameterReader</code> class}
+
+                                // 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.
 class ParameterReader : public Subscriptor
 {
   public:
@@ -159,43 +175,43 @@ class ParameterReader : public Subscriptor
     ParameterHandler &prm;
 };
 
-                               // The constructor stores a reference to 
-                               // the ParameterHandler object that is passed to it:
+                                // The constructor stores a reference to 
+                                // the ParameterHandler object that is passed to it:
 ParameterReader::ParameterReader(ParameterHandler &paramhandler)
                :
                prm(paramhandler)
 {}
 
-                               // @sect4{<code>ParameterReader::declare_parameters</code>}
-
-                               // 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:
+                                // @sect4{<code>ParameterReader::declare_parameters</code>}
+
+                                // 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:
 void ParameterReader::declare_parameters()
 {
-                               // Parameters for mesh and geometry
-                               // include the number of global
-                               // refinement steps that are applied
-                               // to the initial coarse mesh and the
-                               // focal distance $d$ of the
-                               // transducer lens. For the number of
-                               // refinement steps, we allow integer
-                               // values in the range $[0,\infty)$,
-                               // 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:
+                                  // Parameters for mesh and geometry
+                                  // include the number of global
+                                  // refinement steps that are applied
+                                  // to the initial coarse mesh and the
+                                  // focal distance $d$ of the
+                                  // transducer lens. For the number of
+                                  // refinement steps, we allow integer
+                                  // values in the range $[0,\infty)$,
+                                  // 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:
   prm.enter_subsection ("Mesh & geometry parameters");
   {
     prm.declare_entry("Number of refinements", "6",
@@ -210,15 +226,15 @@ void ParameterReader::declare_parameters()
   }
   prm.leave_subsection ();
 
-                               // The next subsection is devoted to
-                               // the physical parameters appearing
-                               // in the equation, which are the
-                               // frequency $\omega$ and wave speed
-                               // $c$. Again, both need to lie in the
-                               // half-open interval $[0,\infty)$
-                               // represented by calling the
-                               // Patterns::Double class with only
-                               // the left end-point as argument:
+                                  // The next subsection is devoted to
+                                  // the physical parameters appearing
+                                  // in the equation, which are the
+                                  // frequency $\omega$ and wave speed
+                                  // $c$. Again, both need to lie in the
+                                  // half-open interval $[0,\infty)$
+                                  // represented by calling the
+                                  // Patterns::Double class with only
+                                  // the left end-point as argument:
   prm.enter_subsection ("Physical constants");
   {
     prm.declare_entry("c", "1.5e5",
@@ -232,50 +248,88 @@ void ParameterReader::declare_parameters()
   prm.leave_subsection ();
 
 
-                               // 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:
+                                  // 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:
   prm.enter_subsection ("Output parameters");
   {
     prm.declare_entry("Output file", "solution",
                      Patterns::Anything(),
                      "Name of the output file (without extension)");
 
-                               // 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. 
+                                    // 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.
     DataOutInterface<1>::declare_parameters (prm);
   }
   prm.leave_subsection ();
 }
 
-                               // @sect4{<code>ParameterReader::read_parameters</code>}
-
-                               // 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:
+                                // @sect4{<code>ParameterReader::read_parameters</code>}
+
+                                // 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:
 void ParameterReader::read_parameters (const std::string parameter_file)
 {
   declare_parameters();
@@ -285,65 +339,98 @@ void ParameterReader::read_parameters (const std::string parameter_file)
 
 
 
-                               // @sect3{The <code>ComputeIntensity</code> class}
-
-                               // As mentioned in the introduction, the quantitiy that we 
-                               // are really after is the spatial distribution of 
-                               // the intensity of the ultrasound wave, which corresponds 
-                               // to $|u|=\sqrt{v^2+w^2}$. Now we could just be content with
-                               // having $v$ and $w$ in our output, and use a suitable
-                               // visualization or postprocessing tool to derive $|u|$ 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 compute_derived_quantities 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 $|u|$ is such a quantity and 
-                               // its computation doesn't even involve any derivatives of $v$ or $w$. 
-
-                               // 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.
-                               // This is what the <code>ComputeIntensity</code> class is about. 
-                               // Notice that all its member functions are implementations of 
-                               // virtual functions defined by the interface class DataPostprocessor.
+                                // @sect3{The <code>ComputeIntensity</code> class}
+
+                                // As mentioned in the introduction,
+                                // the quantitiy that we are really
+                                // after is the spatial distribution
+                                // of the intensity of the ultrasound
+                                // wave, which corresponds to
+                                // $|u|=\sqrt{v^2+w^2}$. Now we could
+                                // just be content with having $v$
+                                // and $w$ in our output, and use a
+                                // suitable visualization or
+                                // postprocessing tool to derive
+                                // $|u|$ 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 compute_derived_quantities
+                                // 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 $|u|$ is such
+                                // a quantity and its computation
+                                // doesn't even involve any
+                                // derivatives of $v$ or $w$.
+
+                                // 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.  This
+                                // is what the
+                                // <code>ComputeIntensity</code>
+                                // class is about.  Notice that all
+                                // its member functions are
+                                // implementations of virtual
+                                // functions defined by the interface
+                                // class DataPostprocessor.
 template <int dim>
 class ComputeIntensity : public DataPostprocessor<dim>
 {
   public:
 
     void compute_derived_quantities_vector (
-                       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 > > &,
-                       std::vector< Vector< double > > &
-                       ) const;
+      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 > > &,
+      std::vector< Vector< double > > &
+    ) const;
 
     std::vector<std::string> get_names () const;
     UpdateFlags              get_needed_update_flags () const;
     unsigned int             n_output_variables () const;
 };
 
-                               // The <code>get_names</code> function returns a vector of strings 
-                               // representing the names we assign to the individual
-                               // quantities that our postprocessor outputs. In our
-                               // case, the postprocessor has only $|u|$ as an output, so we 
-                               // return a vector with a single component named "Intensity":
+                                // The <code>get_names</code>
+                                // function returns a vector of
+                                // strings representing the names we
+                                // assign to the individual
+                                // quantities that our postprocessor
+                                // outputs. In our case, the
+                                // postprocessor has only $|u|$ as an
+                                // output, so we return a vector with
+                                // a single component named
+                                // "Intensity":
 template <int dim>
 std::vector<std::string>
 ComputeIntensity<dim>::get_names() const
@@ -351,19 +438,27 @@ ComputeIntensity<dim>::get_names() const
   return std::vector<std::string> (1, "Intensity");
 }
 
-                               // The next function returns 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 is really needed 
-                               // should be given here, just as we do when we use FEValues objects.
-                               // In our case, only the function values 
-                               // of $v$ and $w$ are needed to compute $|u|$, so we're good 
-                               // with the update_values flag. 
+                                // The next function returns 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 is really needed should
+                                // be given here, just as we do when
+                                // we use FEValues objects.  In our
+                                // case, only the function values of
+                                // $v$ and $w$ are needed to compute
+                                // $|u|$, so we're good with the
+                                // update_values flag.
 template <int dim>
 UpdateFlags
 ComputeIntensity<dim>::get_needed_update_flags () const
@@ -371,11 +466,13 @@ ComputeIntensity<dim>::get_needed_update_flags () const
   return update_values;
 }
 
-                               // To allow the caller to find out how many derived quantities 
-                               // are returned by the postprocessor, the 
-                               // <code>n_output_variables</code> function is used. Since
-                               // we compute only $|u|$, the correct value to return 
-                               // in our case is just 1:
+                                // To allow the caller to find out
+                                // how many derived quantities are
+                                // returned by the postprocessor, the
+                                // <code>n_output_variables</code>
+                                // function is used. Since we compute
+                                // only $|u|$, the correct value to
+                                // return in our case is just 1:
 template <int dim>
 unsigned int
 ComputeIntensity<dim>::n_output_variables () const
@@ -384,61 +481,82 @@ ComputeIntensity<dim>::n_output_variables () const
 }
 
 
-                               // The actual prostprocessing happens in the following function. 
-                               // Its inputs are a vector representing values of the function
-                               // (which is here vector-valued) representing the data vector
-                               // given to DataOut::add_data_vector, evaluated at all quadrature
-                               // points where we generate output,
-                               // and some tensor objects representing derivatives (that we don't 
-                               // use here since $|u|$ is computed from just $v$ and $w$, and for
-                               // which we assign no name to the corresponding function argument). 
-                               // 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. 
+                                // The actual prostprocessing happens
+                                // in the following function.  Its
+                                // inputs are a vector representing
+                                // values of the function (which is
+                                // here vector-valued) representing
+                                // the data vector given to
+                                // DataOut::add_data_vector,
+                                // evaluated at all quadrature points
+                                // where we generate output, and some
+                                // tensor objects representing
+                                // derivatives (that we don't use
+                                // here since $|u|$ is computed from
+                                // just $v$ and $w$, and for which we
+                                // assign no name to the
+                                // corresponding function argument).
+                                // 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.
 template <int dim>
 void
 ComputeIntensity<dim>::compute_derived_quantities_vector (
-                       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*/,
-                       std::vector< Vector< double > >                       & computed_quantities
-                       ) const
+  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*/,
+  std::vector< Vector< double > >                       & computed_quantities
+) const
 {
   Assert(computed_quantities.size() == uh.size(), 
          ExcDimensionMismatch (computed_quantities.size(), uh.size()));
 
-                               // The computation itself is straightforward: We iterate 
-                               // over each entry in the output vector and compute 
-                               // $|u|$ from the corresponding values of $v$ and $w$:
+                                  // The computation itself is
+                                  // straightforward: We iterate over
+                                  // each entry in the output vector
+                                  // and compute $|u|$ from the
+                                  // corresponding values of $v$ and
+                                  // $w$:
   for (unsigned int i=0; i<computed_quantities.size(); i++)
-  {
-    Assert(computed_quantities[i].size() == 1, 
-           ExcDimensionMismatch (computed_quantities[i].size(), 1));
-    Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));
+    {
+      Assert(computed_quantities[i].size() == 1, 
+            ExcDimensionMismatch (computed_quantities[i].size(), 1));
+      Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));
 
-    computed_quantities[i](0) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
-  }
+      computed_quantities[i](0) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
+    }
 }
 
 
-                               // @sect3{The <code>UltrasoundProblem</code> class}
-
-                               // Finally here is the main class of this program. 
-                               // It's member functions are very similar to the previous 
-                               // examples, in particular 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. 
+                                // @sect3{The <code>UltrasoundProblem</code> class}
+
+                                // Finally here is the main class of
+                                // this program.  It's member
+                                // functions are very similar to the
+                                // previous examples, in particular
+                                // 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.
 template <int dim>
 class UltrasoundProblem 
 {
@@ -467,10 +585,14 @@ class UltrasoundProblem
 
 
 
-                               // 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 $v$ and one for $w$: 
+                                // 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 $v$ and
+                                // one for $w$:
 template <int dim>
 UltrasoundProblem<dim>::UltrasoundProblem (ParameterHandler&  param) 
                :
@@ -486,25 +608,31 @@ UltrasoundProblem<dim>::~UltrasoundProblem ()
   dof_handler.clear();
 }
 
-                               // @sect4{<code>UltrasoundProblem::make_grid</code>}
+                                // @sect4{<code>UltrasoundProblem::make_grid</code>}
 
-                               // 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.
+                                // 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.
 template <int dim>
 void UltrasoundProblem<dim>::make_grid ()
 {
-                               // First we generate some logging output 
-                               // and start a timer so we can
-                               // compute execution time when this function is done:
+                                  // First we generate some logging
+                                  // output and start a timer so we
+                                  // can compute execution time when
+                                  // this function is done:
   deallog << "Generating grid... ";
   Timer timer;
   timer.start ();
 
-                               // Then we query the values for the focal distance of the 
-                               // transducer lens and the number of mesh refinement steps 
-                               // from our ParameterHandler object:
+                                  // Then we query the values for the
+                                  // focal distance of the transducer
+                                  // lens and the number of mesh
+                                  // refinement steps from our
+                                  // ParameterHandler object:
   prm.enter_subsection ("Mesh & geometry parameters");
 
   const double         focal_distance = prm.get_double("Focal distance");
@@ -512,37 +640,55 @@ void UltrasoundProblem<dim>::make_grid ()
 
   prm.leave_subsection ();
 
-                               // 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. We 
-                               // compute the radius of this circle in such a way that the 
-                               // segment fits in the interval [0.4,0.6] on the x-axis. 
-                               // 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:
+                                  // 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. We compute
+                                  // the radius of this circle in
+                                  // such a way that the segment fits
+                                  // in the interval [0.4,0.6] on the
+                                  // x-axis.  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:
   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, focal_distance) :
-                                     Point<dim> (0.5, 0.5, focal_distance);
+                      focal_point = (dim == 2) ?
+                                    Point<dim> (0.5, focal_distance) :
+                                    Point<dim> (0.5, 0.5, focal_distance);
 
   const double radius = std::sqrt( (focal_point.distance(transducer) * 
                                    focal_point.distance(transducer)) + 
                                   ((dim==2) ? 0.01 : 0.02));
 
 
-                               // 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 $[0.4,0.6]$ 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 boundary indicator.
+                                  // 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 $[0.4,0.6]$
+                                  // 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 boundary
+                                  // indicator.
   GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
 
   typename Triangulation<dim>::cell_iterator
@@ -556,31 +702,36 @@ void UltrasoundProblem<dim>::make_grid ()
 
         cell->face(face)->set_boundary_indicator (1);
 
-                               // For the circle part of the
-                               // transducer lens, a hyper-ball
-                               // object is used (which, of course,
-                               // in 2D just represents a circle),
-                               // with radius and center as computed
-                               // above. By marking this object as
-                               // <code>static</code>, we ensure that
-                               // it lives until the end of the
-                               // program and thereby longer than the
-                               // triangulation object we will
-                               // associated with it. We then assign
-                               // this boundary-object to the part of
-                               // the boundary with boundary
-                               // indicator 1:
+                                  // For the circle part of the
+                                  // transducer lens, a hyper-ball
+                                  // object is used (which, of course,
+                                  // in 2D just represents a circle),
+                                  // with radius and center as computed
+                                  // above. By marking this object as
+                                  // <code>static</code>, we ensure that
+                                  // it lives until the end of the
+                                  // program and thereby longer than the
+                                  // triangulation object we will
+                                  // associated with it. We then assign
+                                  // this boundary-object to the part of
+                                  // the boundary with boundary
+                                  // indicator 1:
   static const HyperBallBoundary<dim> boundary(focal_point, radius);
   triangulation.set_boundary(1, boundary);
 
-                               // Now global refinement is executed. Cells near the transducer 
-                               // location will be automatically refined according to the 
-                               // circle shaped boundary of the transducer lens:
+                                  // Now global refinement is
+                                  // executed. Cells near the
+                                  // transducer location will be
+                                  // automatically refined according
+                                  // to the circle shaped boundary of
+                                  // the transducer lens:
   triangulation.refine_global (n_refinements);
 
-                               // 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:
+                                  // 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:
   timer.stop ();
   deallog << "done (" 
          << timer()
@@ -593,13 +744,16 @@ void UltrasoundProblem<dim>::make_grid ()
 } 
 
 
-                               // @sect4{<code>UltrasoundProblem::setup_system</code>}
-                               //
-                               // 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:
+                                // @sect4{<code>UltrasoundProblem::setup_system</code>}
+                                //
+                                // 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:
 template <int dim>
 void UltrasoundProblem<dim>::setup_system ()
 {
@@ -632,9 +786,10 @@ void UltrasoundProblem<dim>::setup_system ()
 }
 
 
-                               // @sect4{<code>UltrasoundProblem::assemble_system</code>}
-                               // As before, this function takes care of assembling the 
-                               // system matrix and right hand side vector:
+                                // @sect4{<code>UltrasoundProblem::assemble_system</code>}
+                                // As before, this function takes
+                                // care of assembling the system
+                                // matrix and right hand side vector:
 template <int dim>
 void UltrasoundProblem<dim>::assemble_system () 
 {
@@ -642,10 +797,12 @@ void UltrasoundProblem<dim>::assemble_system ()
   Timer timer;
   timer.start ();
 
-                               // 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.
+                                  // 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.
 
   prm.enter_subsection ("Physical constants");
 
@@ -654,10 +811,15 @@ void UltrasoundProblem<dim>::assemble_system ()
 
   prm.leave_subsection ();
 
-                               // As usual, for computing integrals ordinary Gauss quadrature
-                               // rule is used. Since our bilinear form involves boundary integrals
-                               // on $\Gamma_2$, we also need a quadrature rule for surface
-                               // integration on the faces, which are $dim-1$ dimensional:
+                                  // As usual, for computing
+                                  // integrals ordinary Gauss
+                                  // quadrature rule is used. Since
+                                  // our bilinear form involves
+                                  // boundary integrals on
+                                  // $\Gamma_2$, we also need a
+                                  // quadrature rule for surface
+                                  // integration on the faces, which
+                                  // are $dim-1$ dimensional:
   QGauss<dim>    quadrature_formula(2);
   QGauss<dim-1>  face_quadrature_formula(2);
 
@@ -665,12 +827,19 @@ void UltrasoundProblem<dim>::assemble_system ()
                     n_face_q_points  = face_quadrature_formula.n_quadrature_points,
                     dofs_per_cell    = fe.dofs_per_cell;
 
-                               // The FEValues objects will evaluate the shape functions for us. 
-                               // For the part of the bilinear form that involves integration on 
-                               // $\Omega$, 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. 
+                                  // The FEValues objects will
+                                  // evaluate the shape functions for
+                                  // us.  For the part of the
+                                  // bilinear form that involves
+                                  // integration on $\Omega$, 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.
   FEValues<dim>  fe_values (fe, quadrature_formula, 
                            update_values | update_gradients |
                            update_JxW_values);
@@ -678,10 +847,14 @@ void UltrasoundProblem<dim>::assemble_system ()
   FEFaceValues<dim> fe_face_values (fe, face_quadrature_formula, 
                                    update_values | update_JxW_values);
 
-                               // 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. 
+                                  // 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.
   FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
   std::vector<unsigned int> local_dof_indices (dofs_per_cell);
 
@@ -690,158 +863,360 @@ void UltrasoundProblem<dim>::assemble_system ()
     endc = dof_handler.end();
 
   for (; cell!=endc; ++cell)
-  {
-
-                               // 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: 
-    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)
-      {
-
-                               // 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 $\phi$'s or $\psi$'s
-                               // from the introduction, and similarly
-                               // the corresponding degrees of freedom can be attributed to either $\alpha$ or $\beta$. 
-                               // 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 $i$ and $j$ 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 step-8 or step-22, 
-                               // where these topics are explained in depth. 
-        if (fe.system_to_component_index(i).first == 
-            fe.system_to_component_index(j).first)
-        {
-
-                               // If both DoFs $i$ and $j$ belong to same component, i.e. their shape functions are 
-                               // both $\phi$'s or both $\psi$'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 
-                               // $\phi$ or $\psi$ 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. 
-
-          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));
-
-                               // 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. 
-        }
-      }
-    }
 
-
-                               // 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 $\Gamma_2$, the
-                               // part of the boundary where we have
-                               // absorbing boundary conditions:
-    for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
-      if (cell->face(face)->at_boundary() &&
-         (cell->face(face)->boundary_indicator() == 0) )
-        {
-
-
-                               // 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:
-          fe_face_values.reinit (cell, face);
-
-
-                               // 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:
-          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))
-                               // 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 $\psi$ and which one 
-                               // is a $\phi$, 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 $i$ and $j$ belong to different components, it suffices here
-                               // to test for one of them to which component it belongs.
-                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);
-        }
-
-                               // 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:
-        cell->get_dof_indices (local_dof_indices);
-
-
-                               // and then add the entries to the system matrix one by one:
-    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));
-  }
+                                      // 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:
+      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)
+           {
+
+                                              // 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
+                                              // $\phi$'s or $\psi$'s
+                                              // from the
+                                              // introduction, and
+                                              // similarly the
+                                              // corresponding
+                                              // degrees of freedom
+                                              // can be attributed to
+                                              // either $\alpha$ or
+                                              // $\beta$.  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 $i$ and
+                                              // $j$ 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 step-8 or
+                                              // step-22, where these
+                                              // topics are explained
+                                              // in depth.
+             if (fe.system_to_component_index(i).first == 
+                 fe.system_to_component_index(j).first)
+               {
+
+                                                  // If both DoFs $i$
+                                                  // and $j$ belong
+                                                  // to same
+                                                  // component,
+                                                  // i.e. their shape
+                                                  // functions are
+                                                  // both $\phi$'s or
+                                                  // both $\psi$'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 $\phi$ or
+                                                  // $\psi$ 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.
+
+                 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));
+
+                                                  // 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.
+               }
+           }
+       }
+
+
+                                      // 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 $\Gamma_2$, the
+                                      // part of the boundary where we have
+                                      // absorbing boundary conditions:
+      for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
+       if (cell->face(face)->at_boundary() &&
+           (cell->face(face)->boundary_indicator() == 0) )
+         {
+
+
+                                            // 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:
+           fe_face_values.reinit (cell, face);
+
+
+                                            // 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:
+           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))
+                                                  // 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 $\psi$ and
+                                                  // which one is a
+                                                  // $\phi$, 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 $i$ and
+                                                  // $j$ belong to
+                                                  // different
+                                                  // components, it
+                                                  // suffices here to
+                                                  // test for one of
+                                                  // them to which
+                                                  // component it
+                                                  // belongs.
+                 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);
+         }
+
+                                      // 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...
+      cell->get_dof_indices (local_dof_indices);
+
+
+                                      // ...and then add the entries to
+                                      // the system matrix one by
+                                      // one:
+      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));
+    }
 
 
-                               // The only thing left are the Dirichlet boundary values on 
-                               // $\Gamma_1$, which is characterized by the boundary
-                               // indicator 1. The Dirichlet values are provided by 
-                               // the <code>DirichletBoundaryValues</code> class we defined above:
+                                  // The only thing left are the
+                                  // Dirichlet boundary values on
+                                  // $\Gamma_1$, which is
+                                  // characterized by the boundary
+                                  // indicator 1. The Dirichlet
+                                  // values are provided by the
+                                  // <code>DirichletBoundaryValues</code>
+                                  // class we defined above:
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            1,
@@ -862,24 +1237,37 @@ void UltrasoundProblem<dim>::assemble_system ()
 
 
 
-                               // @sect4{<code>UltrasoundProblem::solve</code>}
+                                // @sect4{<code>UltrasoundProblem::solve</code>}
                                
-                               // 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, which 
-                               // can be achieved by providing the <code> --with-umfpack</code>
-                               // switch to the configure script prior to compilation of the library. 
+                                // 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, which can be
+                                // achieved by providing the <code>
+                                // --with-umfpack</code> switch to
+                                // the configure script prior to
+                                // compilation of the library.
 template <int dim>
 void UltrasoundProblem<dim>::solve ()
 {
@@ -887,16 +1275,29 @@ void UltrasoundProblem<dim>::solve ()
   Timer timer;
   timer.start ();
 
-                               // 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. 
+                                  // 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.
   SparseDirectUMFPACK  A_direct;
   A_direct.initialize(system_matrix);
 
-                               // 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:
+                                  // 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:
   A_direct.vmult (solution, system_rhs);
 
   timer.stop ();
@@ -908,15 +1309,21 @@ void UltrasoundProblem<dim>::solve ()
 
 
 
-                               // @sect4{<code>UltrasoundProblem::output_results</code>}
+                                // @sect4{<code>UltrasoundProblem::output_results</code>}
 
-                               // Here we output our solution $v$ and $w$ as well as the 
-                               // derived quantity $|u|$ in the
-                               // format specified in the parameter file. Most of the 
-                               // work for deriving $|u|$ from $v$ and $w$ 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. 
+                                // Here we output our solution $v$
+                                // and $w$ as well as the derived
+                                // quantity $|u|$ in the format
+                                // specified in the parameter
+                                // file. Most of the work for
+                                // deriving $|u|$ from $v$ and $w$
+                                // 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.
 template <int dim>
 void UltrasoundProblem<dim>::output_results () const
 {
@@ -924,20 +1331,27 @@ void UltrasoundProblem<dim>::output_results () const
   Timer timer;
   timer.start ();
 
-                               // Define objects of our <code>ComputeIntensity</code> class and a DataOut
-                               // object:
+                                  // Define objects of our
+                                  // <code>ComputeIntensity</code>
+                                  // class and a DataOut object:
   ComputeIntensity<dim> intensities;
   DataOut<dim> data_out;
 
   data_out.attach_dof_handler (dof_handler);
 
-                               // 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. 
+                                  // 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.
   prm.enter_subsection("Output parameters");
 
   const std::string output_file    = prm.get("Output file"),
@@ -946,39 +1360,51 @@ void UltrasoundProblem<dim>::output_results () const
 
   prm.leave_subsection ();
 
-                               // Since the ParameterHandler provides the output format 
-                               // parameter as a string, we need to convert it to 
-                               // a format flag that can be understood by the DataOut object. 
-                               // The following function takes care of this: 
+                                  // Since the ParameterHandler
+                                  // provides the output format
+                                  // parameter as a string, we need
+                                  // to convert it to a format flag
+                                  // that can be understood by the
+                                  // DataOut object.  The following
+                                  // function takes care of this:
   DataOutBase::OutputFormat  format = DataOutBase::parse_output_format(output_format);
 
-                               // Now we put together the filename from the base name provided 
-                               // by the ParameterHandler and the suffix which is derived 
-                               // from the format by the DataOutBase::default_suffix function:
+                                  // Now we put together the filename
+                                  // from the base name provided by
+                                  // the ParameterHandler and the
+                                  // suffix which is derived from the
+                                  // format by the
+                                  // DataOutBase::default_suffix
+                                  // function:
   const std::string filename = output_file +
                               DataOutBase::default_suffix(format);
 
   std::ofstream output (filename.c_str());
 
-                               // The solution vectors $v$ and $w$ are added to the DataOut
-                               // object in the usual way:
+                                  // The solution vectors $v$ and $w$
+                                  // are added to the DataOut object
+                                  // in the usual way:
   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);
 
-                               // 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 $|u|$ to the output data:
+                                  // 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 $|u|$ to
+                                  // the output data:
   data_out.add_data_vector (solution, intensities);
 
-                               // 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:
+                                  // 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:
   data_out.build_patches ();
   data_out.write (output, format);
 
@@ -991,8 +1417,9 @@ void UltrasoundProblem<dim>::output_results () const
 
 
 
-                               // @sect4{<code>UltrasoundProblem::run</code>}
-                               // Here we simply execute our functions one after the other:
+                                // @sect4{<code>UltrasoundProblem::run</code>}
+                                // Here we simply execute our
+                                // functions one after the other:
 template <int dim>
 void UltrasoundProblem<dim>::run () 
 {
@@ -1004,55 +1431,55 @@ void UltrasoundProblem<dim>::run ()
 }
 
 
-                               // @sect4{The <code>main</code> function}
-
-                               // 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>step-29.prm</code>. The
-                               // values so read are then handed over
-                               // to an instance of the
-                               // UltrasoundProblem class:
+                                // @sect4{The <code>main</code> function}
+
+                                // 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>step-29.prm</code>. The
+                                // values so read are then handed over
+                                // to an instance of the
+                                // UltrasoundProblem class:
 int main () 
 {
   try
-  {
-    ParameterHandler  prm;
-    ParameterReader   param(prm);
-    param.read_parameters("step-29.prm");
+    {
+      ParameterHandler  prm;
+      ParameterReader   param(prm);
+      param.read_parameters("step-29.prm");
 
-    UltrasoundProblem<2>  ultrasound_problem (prm);
-    ultrasound_problem.run ();
-  }
+      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;
-  }
+    {
+      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;
-}
+}
\ No newline at end of file

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.