From: allmaras Date: Thu, 25 Oct 2007 22:57:24 +0000 (+0000) Subject: Updated step-29 with commented source code and results X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=144963aeb3608e42c4db382700af83541b3c957e;p=dealii-svn.git Updated step-29 with commented source code and results git-svn-id: https://svn.dealii.org/trunk@15378 0785d39b-7218-0410-832d-ea1e28bc413d --- diff --git a/deal.II/examples/step-29/doc/intro.dox b/deal.II/examples/step-29/doc/intro.dox index 0360423334..f2fc0d7231 100644 --- a/deal.II/examples/step-29/doc/intro.dox +++ b/deal.II/examples/step-29/doc/intro.dox @@ -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 equations. This short example demonstrates how this can be implemented in deal.II by using an FE_system object to stack two finite element -fields representing real and imaginary parts. We will also revisit the +fields representing real and imaginary parts. We also revisit the ParameterHandler class introduced in @ref step_19 "step-19", which provides a convenient way for reading parameters from a configuration file at runtime without the need to recompile the program code. @@ -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 frequency ${\omega}>0$ and constant amplitude (that we chose to be 1 here): @f[ -U(x,t) = \sin{\omega t}, \qquad x\in \Gamma_1 +U(x,t) = \cos{\omega t}, \qquad x\in \Gamma_1 @f] If there are no other (interior or boundary) sources, and since the only @@ -147,39 +147,38 @@ c\omega \langle \psi_i,\phi_j\rangle_{\mathrm{L}^2(\Gamma_2)} & -\omega^2 \langl \end{array} \right) @f] -(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.) -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. +(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.) +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. Of course, there is no necessity to choose the spaces $V_h$ and $W_h$ to be -the same. However, we expect the real and imaginary parts of the solution to -have similar properties and will therefore indeed choose $V_h=W_h$. We choose -the current notation using different symbols merely to be able to distinguish -between test functions for $v$ and $w$, as this distinction plays an important -role in the implementation. +the same. However, we expect real and imaginary part of the solution to +have similar properties and will therefore indeed take $V_h=W_h$ in the +implementation, and also use the same basis functions $\phi_i = \psi_i$ for +both spaces. The reason for the notation using different symbols is just that +it allows us to distinguish between shape functions for $v$ and $w$, as this +distinction plays an important role in the implementation.

The test case

-For the computations, we will consider wave propagation in a rectangular -area. Sound is generated by a transducer that we model using a part of the -boundary shaped like a segment of the circle with center at $(0.5, d)$ and a -radius slightly greater than $d$; this shape leads to a focusing of the sound -wave to the center of the circle. Consequently, we take $\Omega$ to be the -unit square $[0,1]^2$ with the piece $[0.4,0.6]\times\{0\}$ removed and -replaced by the transducer lens (see the Results section below for a graphical -depiction of this domain). We let $d$ vary to change the "focus" of the lens -and see how that affects the spatial distribution of the amplitude of $u$, -i.e. we will consider how well focussed $|u|=\sqrt{v^2+w^2}$ is. +For the computations, we will consider wave propagation in the unit square, +with ultrasound generated by a transducer lens that is shaped like a segment +of the circle with center at $(0.5, d)$ and a +radius slightly greater than $d$; this shape should lead to a focusing of the sound +wave at the center of the circle. Varying $d$ changes the "focus" of the lens +and affects the spatial distribution of the intensity of $u$, where our main +concern is how well $|u|=\sqrt{v^2+w^2}$ is.focussed. In the program below, we will implement the complex-valued Helmholtz equations using the formulation with split real and imaginary parts. We will also discuss how to generate a domain that looks like a square with a slight bulge simulating the transducer (in the UltrasoundProblem::make_grid() function), and how to -generate graphical output that not only contains the solution components $v,w$ -but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in -UltrasoundProblem::output_results()). Finally, we will -explain how to use the ParameterHandler class to easily handle the situation -where we want to prescribe the focal distance $d$, wave speed $c$, frequency -$\omega$, and a number of other parameters, in an input file that is read at -run-time, rather than in the source code where we would have to re-compile the -file every time we want to change parameters. +generate graphical output that not only contains the solution components $v$ and +$w$, but also the magnitude $\sqrt{v^2+w^2}$ directly in the output file (in +UltrasoundProblem::output_results()). Finally, we use the +ParameterHandler class to easily read parameters like the focal distance $d$, +wave speed $c$, frequency $\omega$, and a number of other parameters from an +input file at run-time, rather than fixing those parameters in the source code +where we would have to re-compile every time we want to change parameters. diff --git a/deal.II/examples/step-29/doc/results.dox b/deal.II/examples/step-29/doc/results.dox index 747bdbbb3d..ef1f7b8e56 100644 --- a/deal.II/examples/step-29/doc/results.dox +++ b/deal.II/examples/step-29/doc/results.dox @@ -1,3 +1,58 @@

Results

+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 index 0000000000..c6d9609eaf 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 index 0000000000..2b7dc06334 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 index 0000000000..ad29d9d247 Binary files /dev/null and b/deal.II/examples/step-29/doc/step-29.w.png differ diff --git a/deal.II/examples/step-29/step-29.cc b/deal.II/examples/step-29/step-29.cc index 7b9547486e..a9a2bad2f1 100644 --- a/deal.II/examples/step-29/step-29.cc +++ b/deal.II/examples/step-29/step-29.cc @@ -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 @@ -30,20 +40,53 @@ #include - // 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 + // 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 --with-umfpack + // switch when configuring the library: #include + + // 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 + // 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 + + // 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 ctime 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 DirichletBoundaryValues 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 class DirichletBoundaryValues : public Function { @@ -58,6 +101,11 @@ class DirichletBoundaryValues : public Function }; + // 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 inline void DirichletBoundaryValues::vector_value (const Point &/*p*/, @@ -81,36 +129,47 @@ void DirichletBoundaryValues::vector_value_list (const std::vector::vector_value (points[p], value_list[p]); } + // @sect3{The ParameterReader class} + // The next classis responsible for preparing the + // ParameterHandler object and reading parameters from + // an input file. + // It includes a function declare_parameters + // that declares all the necessary parameters + // and a read_parameters + // 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 ¶mhandler) : prm(paramhandler) {} + // @sect4{ParameterReader::declare_parameters} -void ParameterReader::read_parameters() -{ - declare_parameters(); - - const std::string parameter_file = "step-29.prm"; - prm.read_input (parameter_file); -} - - + // The declare_parameters 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 FormatFlags::declare_parameters + // function, which declares all the parameters specific to that format in + // an own subsection. The following call of + // DataOutInterface<1>::declare_parameters executes + // declare_parameters 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{ParameterReader::read_parameters} + + // 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 ComputeIntensity 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 ComputeIntensity class is about. + // Notice that all its member functions are implementations of + // virtual functions defined by the interface class DataPostprocessor. template -class Postprocessor : public DataPostprocessor +class ComputeIntensity : public DataPostprocessor { public: @@ -158,7 +292,7 @@ class Postprocessor : public DataPostprocessor 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 unsigned int n_output_variables () const; }; - + // The get_names 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 std::vector -Postprocessor::get_names() const +ComputeIntensity::get_names() const { std::vector 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 UpdateFlags -Postprocessor::get_needed_update_flags () const +ComputeIntensity::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 + // n_output_variables function is used. Since + // we compute only $|u|$, the correct value to return + // in our case is just 1: template unsigned int -Postprocessor::n_output_variables () const +ComputeIntensity::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 + // computed_quantities vector. + // Remember that this function may only use data for which the + // respective update flag is specified by + // get_needed_update_flags. For example, we may + // not use the derivatives here, + // since our implementation of get_needed_update_flags + // requests that only function values are provided. template void -Postprocessor::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::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; iUltrasoundProblem 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 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 system_matrix; Vector solution, system_rhs; - - const double c, omega; }; -template -double -UltrasoundProblem::get_omega (ParameterHandler &prm) -{ - prm.enter_subsection ("Physical constants"); - double omega_tmp = prm.get_double("omega"); - prm.leave_subsection (); - - return omega_tmp; -} - - -template -double -UltrasoundProblem::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 UltrasoundProblem::UltrasoundProblem (ParameterHandler& param) : prm(param), dof_handler(triangulation), - fe(FE_Q(1), 2), - c(get_c(prm)), - omega(get_omega(prm)) + fe(FE_Q(1), 2) {} @@ -297,10 +445,24 @@ UltrasoundProblem::~UltrasoundProblem () dof_handler.clear(); } + // @sect4{UltrasoundProblem::make_grid} + // 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 void UltrasoundProblem::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::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 transducer = (dim == 2) ? Point (0.5, 0.0) : Point (0.5, 0.5, 0.0), @@ -321,6 +491,15 @@ void UltrasoundProblem::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::cell_iterator cell = triangulation.begin (), endc = triangulation.end(); @@ -332,27 +511,57 @@ void UltrasoundProblem::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 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 clock() function + // is implementation depended, and also the clock_t 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{UltrasoundProblem::setup_system} + // Initialization of the system matrix, sparsity patterns + // and vectors are the same as in previous examples + // and therefore do not need further comment: template void UltrasoundProblem::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::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{UltrasoundProblem::assemble_system} + // As before, this function takes care of assembling the + // system matrix and right hand side vector: template void UltrasoundProblem::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 quadrature_formula(2); QGauss face_quadrature_formula(2); @@ -380,10 +618,12 @@ void UltrasoundProblem::assemble_system () n_face_q_points = face_quadrature_formula.n_quadrature_points, dofs_per_cell = fe.dofs_per_cell; - FullMatrix cell_matrix (dofs_per_cell, dofs_per_cell); - - std::vector 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 fe_values (fe, quadrature_formula, update_values | update_gradients | update_JxW_values); @@ -391,12 +631,23 @@ void UltrasoundProblem::assemble_system () FEFaceValues 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 cell_matrix (dofs_per_cell, dofs_per_cell); + std::vector local_dof_indices (dofs_per_cell); + typename DoFHandler::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::assemble_system () { for (unsigned int j=0; jat_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::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::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::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::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 DirichletBoundaryValues class we defined above: std::map boundary_values; VectorTools::interpolate_boundary_values (dof_handler, 1, @@ -462,49 +796,142 @@ void UltrasoundProblem::assemble_system () system_matrix, solution, system_rhs); + + clock_t end = clock(); + deallog << "done (" + << (end - start) / tps + << "s)" + << std::endl; } + + // @sect4{UltrasoundProblem::solve} + template void UltrasoundProblem::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 --with-umfpack + // switch to the configure script prior to compilation of the library. SparseDirectUMFPACK A_direct; + // The initialize 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 A_direct 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{UltrasoundProblem::output_results} + + // 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 ComputeIntensity class, + // so that the output routine is rather straightforward and very similar + // to what is done in the previous tutorials. template void UltrasoundProblem::output_results () const { - Postprocessor pproc; + deallog << "Generating output... "; + clock_t start = clock(); + + // Define objects of our ComputeIntensity class and a DataOut + // object: + ComputeIntensity intensities; DataOut 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 + // ParameterReader::declare_parameters. 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 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 add_data_vector again, + // but this with our ComputeIntensity 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{UltrasoundProblem::run} + // Here we simply execute our functions one after the other: template void UltrasoundProblem::run () { @@ -516,17 +943,27 @@ void UltrasoundProblem::run () } + // @sect4{The main function} + + // Finally the main 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 ParameterReader objects, + // and let the latter read in the parameter values from + // a textfile called step-29.prm: + 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 ultrasound_problem (prm); ultrasound_problem.run (); } @@ -553,6 +990,5 @@ int main () << std::endl; return 1; } - return 0; } diff --git a/deal.II/examples/step-29/step-29.prm b/deal.II/examples/step-29/step-29.prm index 9f48aef3ca..13b589c477 100644 --- a/deal.II/examples/step-29/step-29.prm +++ b/deal.II/examples/step-29/step-29.prm @@ -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