]> https://gitweb.dealii.org/ - dealii-svn.git/commitdiff
Updated step-29 with commented source code and results
authorallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Oct 2007 22:57:24 +0000 (22:57 +0000)
committerallmaras <allmaras@0785d39b-7218-0410-832d-ea1e28bc413d>
Thu, 25 Oct 2007 22:57:24 +0000 (22:57 +0000)
git-svn-id: https://svn.dealii.org/trunk@15378 0785d39b-7218-0410-832d-ea1e28bc413d

deal.II/examples/step-29/doc/intro.dox
deal.II/examples/step-29/doc/results.dox
deal.II/examples/step-29/doc/step-29.intensity.png [new file with mode: 0644]
deal.II/examples/step-29/doc/step-29.v.png [new file with mode: 0644]
deal.II/examples/step-29/doc/step-29.w.png [new file with mode: 0644]
deal.II/examples/step-29/step-29.cc
deal.II/examples/step-29/step-29.prm

index 0360423334055bfa6ceb2ac642ce0e14ac4a580c..f2fc0d7231411f8836eb011da79c8222135b557c 100644 (file)
@@ -16,7 +16,7 @@ element fields for discretizing each one of them. Basically this amounts to
 viewing a single complex valued equation as a system of two real valued\r
 equations. This short example demonstrates how this can be implemented in\r
 deal.II by using an <code>FE_system</code> object to stack two finite element\r
-fields representing real and imaginary parts. We will also revisit the\r
+fields representing real and imaginary parts. We also revisit the\r
 ParameterHandler class introduced in @ref step_19 "step-19", which provides a\r
 convenient way for reading parameters from a configuration file at runtime\r
 without the need to recompile the program code.  \r
@@ -47,7 +47,7 @@ choose boundary conditions on $\Gamma_2$ in such a way that they imitate a
 larger domain). On $\Gamma_1$, the transducer generates a wave of constant\r
 frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here):  \r
 @f[\r
-U(x,t) = \sin{\omega t}, \qquad x\in \Gamma_1\r
+U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1\r
 @f]\r
 \r
 If there are no other (interior or boundary) sources, and since the only\r
@@ -147,39 +147,38 @@ c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)} & -\omega^2 \langl
 \end{array}\r
 \right)\r
 @f]\r
-(One should not be fooled by the right hand side being zero here, that is because we haven't included the Dirichlet boundary data yet.) \r
-Because of the alternating sign in the off-diagonal blocks, we can already see that this system is non-symmetric, in fact it is even indefinite.\r
+(One should not be fooled by the right hand side being zero here, that is \r
+because we haven't included the Dirichlet boundary data yet.) \r
+Because of the alternating sign in the off-diagonal blocks, we can already \r
+see that this system is non-symmetric, in fact it is even indefinite.\r
 Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be\r
-the same. However, we expect the real and imaginary parts of the solution to\r
-have similar properties and will therefore indeed choose $V_h=W_h$. We choose\r
-the current notation using different symbols merely to be able to distinguish\r
-between test functions for $v$ and $w$, as this distinction plays an important\r
-role in the implementation.  \r
+the same. However, we expect real and imaginary part of the solution to\r
+have similar properties and will therefore indeed take $V_h=W_h$ in the \r
+implementation, and also use the same basis functions $\phi_i = \psi_i$ for \r
+both spaces. The reason for the notation using different symbols is just that \r
+it allows us to distinguish between shape functions for $v$ and $w$, as this \r
+distinction plays an important role in the implementation.  \r
 \r
 \r
 <h3>The test case</h3>\r
 \r
-For the computations, we will consider wave propagation in a rectangular\r
-area. Sound is generated by a transducer that we model using a part of the\r
-boundary shaped like a segment of the circle with center at $(0.5, d)$ and a\r
-radius slightly greater than $d$; this shape leads to a focusing of the sound\r
-wave to the center of the circle. Consequently, we take $\Omega$ to be the\r
-unit square $[0,1]^2$ with the piece $[0.4,0.6]\times\{0\}$ removed and\r
-replaced by the transducer lens (see the Results section below for a graphical\r
-depiction of this domain). We let $d$ vary to change the "focus" of the lens\r
-and see how that affects the spatial distribution of the amplitude of $u$,\r
-i.e. we will consider how well focussed $|u|=\sqrt{v^2+w^2}$ is.\r
+For the computations, we will consider wave propagation in the unit square, \r
+with ultrasound generated by a transducer lens that is shaped like a segment \r
+of the circle with center at $(0.5, d)$ and a \r
+radius slightly greater than $d$; this shape should lead to a focusing of the sound\r
+wave at the center of the circle. Varying $d$ changes the "focus" of the lens\r
+and affects the spatial distribution of the intensity of $u$, where our main\r
+concern is how well $|u|=\sqrt{v^2+w^2}$ is.focussed. \r
 \r
 In the program below, we will implement the complex-valued Helmholtz equations\r
 using the formulation with split real and imaginary parts. We will also\r
 discuss how to generate a domain that looks like a square with a slight bulge\r
 simulating the transducer (in the\r
 <code>UltrasoundProblem<dim>::make_grid()</code> function), and how to\r
-generate graphical output that not only contains the solution components $v,w$\r
-but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in\r
-<code>UltrasoundProblem<dim>::output_results()</code>). Finally, we will\r
-explain how to use the ParameterHandler class to easily handle the situation\r
-where we want to prescribe the focal distance $d$, wave speed $c$, frequency\r
-$\omega$, and a number of other parameters, in an input file that is read at\r
-run-time, rather than in the source code where we would have to re-compile the\r
-file every time we want to change parameters.\r
+generate graphical output that not only contains the solution components $v$ and \r
+$w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in\r
+<code>UltrasoundProblem<dim>::output_results()</code>). Finally, we use the \r
+ParameterHandler class to easily read parameters like the focal distance $d$, \r
+wave speed $c$, frequency $\omega$, and a number of other parameters from an \r
+input file at run-time, rather than fixing those parameters in the source code \r
+where we would have to re-compile every time we want to change parameters.\r
index 747bdbbb3ddc12b135f3b03463007a96ca80561b..ef1f7b8e5685f7c9dd658fbefcd89132bf56ba82 100644 (file)
@@ -1,3 +1,58 @@
 <a name="Results"></a>
 <h1>Results</h1>
 
+The parameter file that comes with the tutorial code sets
+d=0.3, which amounts to a focus of the transducer lens 
+at x=0.5, y=0.3. The coarse mesh is refined 5 times, 
+resulting in 160x160 cells, and the output is written in gmv 
+format.
+
+Here's the console output of the program in debug mode: 
+
+@code 
+examples/step-29> make run
+============================ Running step-29
+DEAL::Generating grid... done (1.34000s)
+DEAL::  Number of active cells:  25600
+DEAL::Setting up system... done (1.17000s)
+DEAL::  Number of degrees of freedom: 51842
+DEAL::Assembling system matrix... done (7.52000s)
+DEAL::Solving linear system... done (4.43000s)
+DEAL::Generating output... done (5.94000s)
+@endcode
+
+(Of course, execution times will differ if you run the program 
+locally.) The fact that most of the time is spent on assembling 
+the system matrix and generating output is due to some extra checks 
+done in debug mode, in optimized mode these parts of the program  
+run much faster: 
+
+@code
+examples/step-29> make run
+============================ Running step-29
+DEAL::Generating grid... done (0.0500000s)
+DEAL::  Number of active cells:  25600
+DEAL::Setting up system... done (0.170000s)
+DEAL::  Number of degrees of freedom: 51842
+DEAL::Assembling system matrix... done (0.370000s)
+DEAL::Solving linear system... done (3.65000s)
+DEAL::Generating output... done (1.23000s)
+@endcode
+
+The graphical output of the program looks as follows: 
+
+
+@image html step-29.v.png "v = Re(u)" width=5cm
+
+@image html step-29.w.png "w = Im(u)" width=5cm
+
+@image html step-29.intensity.png "|u|" width=5cm
+
+The first two picturse show the real and imaginary parts of 
+$u$, whereas the last shows the intensity $|u|$. One can clearly 
+see that the intensity is focussed around the focal point of the 
+lens (0.5, 0.3), and that the focus 
+is rather sharp in x-direction but more blurred in y-direction, which is a 
+consequence of the geometry and the wave nature of the problem. 
+
+
diff --git a/deal.II/examples/step-29/doc/step-29.intensity.png b/deal.II/examples/step-29/doc/step-29.intensity.png
new file mode 100644 (file)
index 0000000..c6d9609
Binary files /dev/null and b/deal.II/examples/step-29/doc/step-29.intensity.png differ
diff --git a/deal.II/examples/step-29/doc/step-29.v.png b/deal.II/examples/step-29/doc/step-29.v.png
new file mode 100644 (file)
index 0000000..2b7dc06
Binary files /dev/null and b/deal.II/examples/step-29/doc/step-29.v.png differ
diff --git a/deal.II/examples/step-29/doc/step-29.w.png b/deal.II/examples/step-29/doc/step-29.w.png
new file mode 100644 (file)
index 0000000..ad29d9d
Binary files /dev/null and b/deal.II/examples/step-29/doc/step-29.w.png differ
index 7b9547486e2491df86a3448ca0943d9907cda4c7..a9a2bad2f14115425eab0c26f34204f787f04946 100644 (file)
@@ -1,3 +1,13 @@
+/*                                                                */
+/*    Copyright (C) 2007 by the deal.II authors and M. Allmaras   */
+/*                                                                */
+/*    This file is subject to QPL and may not be  distributed     */
+/*    without copyright and license information. Please refer     */
+/*    to the file deal.II/doc/license.html for the text and       */
+/*    further information on this license.                        */
+
+
+
                                // @sect3{Include files}
 
                                // The following header files are unchanged 
 
 #include <fstream>
 
-                               // This header file is needed for  
-                               // the ParameterHandler class that we will
-                               // use to read parameters from a 
-                               // configuration file during runtime:
+                               // 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:
 #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:
 #include <fe/fe_system.h>
 
+                               // The last header is from the C++ standard library and provides
+                               // functions that will allow us to measure execution time
+                               // of various parts of our program:
+#include <ctime>
+
+                               // Although we'll follow good deal.ii practice and keep 
+                               // all of the code dimension independent, we will 
+                               // really only consider the 2D problem here: 
 #define DIM 2
 
 using namespace dealii;
 
+                               // The next line provides a shorthand to the CLOCKS_PER_SEC 
+                               // constant, which is defined in the <code>ctime</code> header
+                               // and contains the number of processor ticks in each second. 
+                               // Henceforth we can compute actual time from ticks by simply dividing
+                               // by tps. 
+const double tps = CLOCKS_PER_SEC;
+
 
+                               // @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. 
 template <int dim>
 class DirichletBoundaryValues : public Function<dim>
 {
@@ -58,6 +101,11 @@ class DirichletBoundaryValues : public Function<dim>
 };
 
 
+                               // Since there are two values $v$ and $w$ that need to be prescribed 
+                               // at the boundary, the boundary value function must return a vector
+                               // 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>
 inline
 void DirichletBoundaryValues<dim>::vector_value (const Point<dim> &/*p*/,
@@ -81,36 +129,47 @@ 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 classis 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:
     ParameterReader(ParameterHandler &);
-    void read_parameters();
+    void read_parameters(const std::string);
 
   private:
     void declare_parameters();
     ParameterHandler &prm;
 };
 
-
+                               // 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>}
 
-void ParameterReader::read_parameters()
-{
-  declare_parameters();
-
-  const std::string parameter_file = "step-29.prm";
-  prm.read_input (parameter_file);
-}
-
-
+                               // The <code>declare_parameters</code> function declares all the parameters 
+                               // that our ParameterHandler object will discover in the input file, 
+                               // along with their types, range conditions and the subsections
+                               // they appear in: 
 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 between 1 and 10, 
+                               // and for the focal distance any number greater than zero:
   prm.enter_subsection ("Mesh & geometry parameters");
 
     prm.declare_entry("Number of refinements", "6",
@@ -121,36 +180,111 @@ void ParameterReader::declare_parameters()
      prm.declare_entry("Focal distance", "0.3",
                       Patterns::Double(0),
                       "Distance of the focal point of the lens "
-                      "to the x-axis (or xy-plane in 3D)");
+                      "to the x-axis");
 
   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$:
   prm.enter_subsection ("Physical constants");
 
     prm.declare_entry("c", "1.5e5",
-                     Patterns::Double(),
+                     Patterns::Double(0),
                      "Wave speed");
 
     prm.declare_entry("omega", "5.0e7",
-                     Patterns::Double(),
+                     Patterns::Double(0),
                      "Frequency");
 
   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:
   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. 
+                               // 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 and first initiates declaration of 
+                               // the parameters, and then tries to read them from the input file whose 
+                               // filename is provided by the caller. 
+void ParameterReader::read_parameters(const std::string parameter_file)
+{
+  declare_parameters();
+
+  prm.read_input (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 contaning 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 Postprocessor : public DataPostprocessor<dim>
+class ComputeIntensity : public DataPostprocessor<dim>
 {
   public:
 
@@ -158,7 +292,7 @@ class Postprocessor : public DataPostprocessor<dim>
                        const std::vector< Vector< double > > &, 
                        const std::vector< std::vector< Tensor< 1, dim > > > &, 
                        const std::vector< std::vector< Tensor< 2, dim > > > &, 
-                       const std::vector< Point< dim > >                    &,
+                       const std::vector< Point< dim > > &,
                        std::vector< Vector< double > > &
                        ) const;
 
@@ -167,63 +301,103 @@ class Postprocessor : public DataPostprocessor<dim>
     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 
+                               // need to provide just that one name:
 template <int dim>
 std::vector<std::string>
-Postprocessor<dim>::get_names() const
+ComputeIntensity<dim>::get_names() const
 {
   std::vector<std::string> field_names;
-
-  field_names.push_back("Re_u");
-  field_names.push_back("Im_u");
   field_names.push_back("Intensity");
 
   return field_names;
 }
 
-
+                               // 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). 
+                               // Of course, computation of the derivatives requires additional 
+                               // resources, so only the flags for data that is really needed 
+                               // should be given here. 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
-Postprocessor<dim>::get_needed_update_flags () const
+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:
 template <int dim>
 unsigned int
-Postprocessor<dim>::n_output_variables () const
+ComputeIntensity<dim>::n_output_variables () const
 {
-  return 3;
+  return 1;
 }
 
 
+                               // The actual prostprocessing happens in the following function. 
+                               // Its inputs are a vector representing point values of the function 
+                               // and some tensor objects representing derivatives (that we don't 
+                               // use here since $|u|$ is computed from just $v$ and $w$). 
+                               // 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
-Postprocessor<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
+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
 {
   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$:
   for (unsigned int i=0; i<computed_quantities.size(); i++)
   {
-    Assert(computed_quantities[i].size() == 3
-           ExcDimensionMismatch (computed_quantities[i].size(), 3));
+    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) = uh[i](0);
-    computed_quantities[i](1) = uh[i](1);
-    computed_quantities[i](2) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
+    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 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 
 {
@@ -233,9 +407,6 @@ class UltrasoundProblem
     void run ();
 
   private:
-    static double get_omega (ParameterHandler &);
-    static double get_c (ParameterHandler &);
-
     void make_grid ();
     void setup_system ();
     void assemble_system ();
@@ -251,43 +422,20 @@ class UltrasoundProblem
     SparsityPattern        sparsity_pattern;
     SparseMatrix<double>   system_matrix;      
     Vector<double>         solution, system_rhs;
-
-    const double           c, omega;
 };
 
 
-template <int dim>
-double
-UltrasoundProblem<dim>::get_omega (ParameterHandler &prm)
-{
-  prm.enter_subsection ("Physical constants");
-    double omega_tmp = prm.get_double("omega");
-  prm.leave_subsection ();
-
-  return omega_tmp;
-}
-
-
-template <int dim>
-double
-UltrasoundProblem<dim>::get_c (ParameterHandler &prm)
-{
-  prm.enter_subsection ("Physical constants");
-    double c_tmp = prm.get_double("c");
-  prm.leave_subsection ();
-
-  return c_tmp;
-}
-
 
+                               // 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) 
   :
   prm(param),
   dof_handler(triangulation),
-  fe(FE_Q<dim>(1), 2),
-  c(get_c(prm)),
-  omega(get_omega(prm))
+  fe(FE_Q<dim>(1), 2)
 {}
 
 
@@ -297,10 +445,24 @@ UltrasoundProblem<dim>::~UltrasoundProblem ()
   dof_handler.clear();
 }
 
+                               // @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 
+                               // 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 store the current number of ticks to be able to 
+                               // compute execution time when this function is done:
+  deallog << "Generating grid... ";
+  clock_t start = clock();
+
+                               // 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");
@@ -308,8 +470,16 @@ void UltrasoundProblem<dim>::make_grid ()
 
   prm.leave_subsection ();
 
-  GridGenerator::subdivided_hyper_cube (triangulation, 5, 0, 1);
-
+                               // 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), 
@@ -321,6 +491,15 @@ void UltrasoundProblem<dim>::make_grid ()
                         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. 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
     cell = triangulation.begin (),
     endc = triangulation.end();
@@ -332,27 +511,57 @@ 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. Then we assign this boundary-object
+                               // to the part of the boundary with boundary indicator 1:
   const HyperBallBoundary<dim> boundary(focal_point, radius);
   triangulation.set_boundary(1, boundary);
 
+                               // Now the 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_ref);
 
-  deallog << " Number of active cells:  "
-         << triangulation.n_active_cells()
+                               // The next line releases the triangulation's 
+                               // pointer to the boundary object that we just created, which 
+                               // is necessary since the boundary object will be destructed 
+                               // as we leave this function 
+                               // and we don't want the triangulation to keep a hanging pointer. 
+  triangulation.set_boundary(1);
+
+                               // Lastly, we generate some more logging output. By querying 
+                               // the present number of ticks again and comparing to 
+                               // what we had at the beginning of the function, we can 
+                               // calculate execution time by dividing by tps. 
+                               // Note that the resolution of the <code>clock()</code> function 
+                               // is implementation depended, and also the <code>clock_t</code> values 
+                               // it returns may overflow, so this way of measuring execution 
+                               // time should be taken with a grain of salt as it may not 
+                               // be very accurate and even completely wrong for longer timespans:
+  clock_t end = clock();
+  deallog << "done (" 
+         << (end - start) / tps 
+         << "s)" 
          << std::endl;
 
-  triangulation.set_boundary(1);
+  deallog << "  Number of active cells:  "
+         << triangulation.n_active_cells()
+         << std::endl;
 } 
 
 
+                               // @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:
 template <int dim>
 void UltrasoundProblem<dim>::setup_system ()
 {
-  dof_handler.distribute_dofs (fe);
+  deallog << "Setting up system... ";
+  clock_t start = clock();
 
-  deallog << " Number of degrees of freedom: "
-         << dof_handler.n_dofs()
-         << std::endl;
+  dof_handler.distribute_dofs (fe);
 
   sparsity_pattern.reinit (dof_handler.n_dofs(),
                           dof_handler.n_dofs(),
@@ -364,15 +573,44 @@ void UltrasoundProblem<dim>::setup_system ()
   system_matrix.reinit (sparsity_pattern);
   system_rhs.reinit (dof_handler.n_dofs());
   solution.reinit (dof_handler.n_dofs());
+
+  clock_t end = clock();
+  deallog << "done (" 
+         << (end - start) / tps 
+         << "s)" 
+         << std::endl;
+
+  deallog << "  Number of degrees of freedom: "
+         << dof_handler.n_dofs()
+         << std::endl;
 }
 
 
+                               // @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 () 
 {
-  const double om2 = omega * omega;
-  const double c2 = c * c;
+  deallog << "Assembling system matrix... ";
+  clock_t start = clock();
+
+                               // 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");
+
+  const double omega = prm.get_double("omega"),
+              c     = prm.get_double("c");
+
+  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:
   QGauss<dim>    quadrature_formula(2);
   QGauss<dim-1>  face_quadrature_formula(2);
 
@@ -380,10 +618,12 @@ void UltrasoundProblem<dim>::assemble_system ()
                     n_face_q_points  = face_quadrature_formula.n_quadrature_points,
                     dofs_per_cell    = fe.dofs_per_cell;
 
-  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
-
-  std::vector<unsigned int> local_dof_indices (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. 
   FEValues<dim>  fe_values (fe, quadrature_formula, 
                            update_values | update_gradients |
                            update_JxW_values);
@@ -391,12 +631,23 @@ 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. 
+  FullMatrix<double> cell_matrix (dofs_per_cell, dofs_per_cell);
+  std::vector<unsigned int> local_dof_indices (dofs_per_cell);
+
   typename DoFHandler<dim>::active_cell_iterator
     cell = dof_handler.begin_active(),
     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);
 
@@ -404,29 +655,89 @@ void UltrasoundProblem<dim>::assemble_system ()
     {
       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, 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. 
+                               // But if you look at the form of the system matrix in the introduction, this disctinction 
+                               // 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, 
+                               // 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)) *
-                                 (- om2)
+                                 (- omega * omega)
                                  +
                                  (fe_values.shape_grad(i,q_point) *
                                   fe_values.shape_grad(j,q_point)) *
-                                 c2) *
+                                 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. 
         }
       }
     }
 
+
+                               // For DoFs that belong to different components of the system, i.e. one DoF
+                               // representing a $\phi$ and the other a $\psi$, a contribution is only 
+                               // possible in the off-diagonal blocks of the system matrix. The entries 
+                               // in these blocks consist of a boundary integral on $\Gamma_2$, so we 
+                               // should first check if the current cell is on the boundary at all, since 
+                               // if it is not, its shape functions will certainly not have support on the boundary. 
     if (cell->at_boundary())
+
+                               // If the current cell is at the boundary, we look through its
+                               // faces to identify the ones that lie on $\Gamma_2$:
       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 != 
@@ -434,8 +745,21 @@ void UltrasoundProblem<dim>::assemble_system ()
                  fe.has_support_on_face(i, face) &&
                  fe.has_support_on_face(j, face))
 
+
+                               // These DoFs will then 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, so 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)) * 
+                  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 *
@@ -443,8 +767,13 @@ void UltrasoundProblem<dim>::assemble_system ()
                                      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],
@@ -452,6 +781,11 @@ void UltrasoundProblem<dim>::assemble_system ()
                           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:
   std::map<unsigned int,double> boundary_values;
   VectorTools::interpolate_boundary_values (dof_handler,
                                            1,
@@ -462,49 +796,142 @@ void UltrasoundProblem<dim>::assemble_system ()
                                      system_matrix,
                                      solution,
                                      system_rhs);
+
+  clock_t end = clock();
+  deallog << "done (" 
+         << (end - start) / tps 
+         << "s)" 
+         << std::endl;
 }
 
 
+
+                               // @sect4{<code>UltrasoundProblem::solve</code>}
+                               
 template <int dim>
 void UltrasoundProblem<dim>::solve ()
 {
+  deallog << "Solving linear system... ";
+  clock_t start = clock();
+
+                               // 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. 
   SparseDirectUMFPACK  A_direct;
 
+                               // The <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. 
   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:
   A_direct.vmult(solution,system_rhs);
+
+  clock_t end = clock();
+  deallog << "done (" 
+         << (end - start) / tps 
+         << "s)" 
+         << std::endl;
 }
 
 
+
+                               // @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. 
 template <int dim>
 void UltrasoundProblem<dim>::output_results () const
 {
-  Postprocessor<dim> pproc;
+  deallog << "Generating output... ";
+  clock_t start = clock();
+
+                               // 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: 
   prm.enter_subsection("Output parameters");
 
   const std::string output_file    = prm.get("Output file"),
                    output_format  = prm.get("Output format");
 
+                               // 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. 
   data_out.parse_parameters(prm);
 
   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: 
   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:
   const std::string filename = output_file +
                               DataOutBase::default_suffix(format);
 
   std::ofstream output (filename.c_str());
 
-  data_out.add_data_vector (solution, pproc);
+                               // 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:
+  data_out.add_data_vector (solution, intensities);
+
+                               // The last steps are as before: 
   data_out.build_patches ();
   data_out.write (output, format);
+
+  clock_t end = clock();
+  deallog << "done (" 
+         << (end - start) / tps 
+         << "s)" 
+         << std::endl;
 }
 
 
+
+                               // @sect4{<code>UltrasoundProblem::run</code>}
+                               // Here we simply execute our functions one after the other:
 template <int dim>
 void UltrasoundProblem<dim>::run () 
 {
@@ -516,17 +943,27 @@ void UltrasoundProblem<dim>::run ()
 }
 
 
+                               // @sect4{The <code>main</code> function}
+
+                               // Finally the <code>main</code> function of the program: 
 int main () 
 {
   try
   {
-    ParameterHandler  prm;
+                               // In 1D, the description of the domain
+                               // and the boundary conditions is not very sensible, so
+                               // exclude this case:
+    Assert (DIM > 1, ExcNotImplemented());
 
+                               // Next 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>:
+    ParameterHandler  prm;
     ParameterReader   param(prm);
-    param.read_parameters();
-
-    Assert (DIM > 1, ExcNotImplemented());
+    param.read_parameters("step-29.prm");
 
+                               // Lastly, we instantiate our main class with the ParameterHandler
+                               // object and start the computations:
     UltrasoundProblem<DIM>  ultrasound_problem (prm);
     ultrasound_problem.run ();
   }
@@ -553,6 +990,5 @@ int main ()
              << std::endl;
     return 1;
   }
-
   return 0;
 }
index 9f48aef3caf1bc92f48bc08b1b82d541c375cf52..13b589c477b448f37296e0567ae1ac5c175c912b 100644 (file)
@@ -2,11 +2,11 @@
 # ---------------------
 
 subsection Mesh & geometry parameters
-  # Distance of the focal point of the lens to the x-axis (or xy-plane in 3D)
+  # Distance of the focal point of the lens to the x-axis
   set Focal distance        = 0.3
 
   # Number of global mesh refinement steps applied to initial coarse grid
-  set Number of refinements = 6
+  set Number of refinements = 5
 end
 
 
@@ -15,7 +15,7 @@ subsection Physical constants
   set c     = 1.5e5
 
   # Frequency
-  set omega = 5.0e7
+  set omega = 3.0e7
 end
 
 

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

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.