From 1157fab4cd8f964716075f83600dc9ecb83d32b9 Mon Sep 17 00:00:00 2001 From: David Wells Date: Sat, 4 May 2019 11:42:57 -0400 Subject: [PATCH] Modernize step-23. This PR introduces several small changes: 1. Don't hard-code two quadrature degrees (followup to #7943). 2. Remove unnecessary default constructor definitions. 3. Combine declarations and definitions of some simple functions> 4. Remove unused headers. 5. Improve whitespace (and remove tabs). 6. Add some more citations to WorkStream. 7. Switch to VTU output so that users can remake the movie. --- examples/step-23/doc/results.dox | 10 +- examples/step-23/step-23.cc | 172 ++++++++++++------------------- 2 files changed, 70 insertions(+), 112 deletions(-) diff --git a/examples/step-23/doc/results.dox b/examples/step-23/doc/results.dox index 11dffad4c0..5146a0d09f 100644 --- a/examples/step-23/doc/results.dox +++ b/examples/step-23/doc/results.dox @@ -115,11 +115,11 @@ If you want to explore a bit, try out some of the following things: possible to avoid copying by not eliminating columns in the linear systems, which is implemented by appending a @p false argument to the call: @code - MatrixTools::apply_boundary_values (boundary_values, - matrix_u, - solution_u, - system_rhs, - false); + MatrixTools::apply_boundary_values(boundary_values, + matrix_u, + solution_u, + system_rhs, + false); @endcode
  • deal.II being a library that supports adaptive meshes it would of course be diff --git a/examples/step-23/step-23.cc b/examples/step-23/step-23.cc index 5b81444550..8708099c77 100644 --- a/examples/step-23/step-23.cc +++ b/examples/step-23/step-23.cc @@ -24,10 +24,8 @@ // many of the previous tests: #include #include -#include #include -#include #include #include #include @@ -36,15 +34,11 @@ #include #include -#include -#include #include -#include #include #include -#include #include @@ -78,7 +72,7 @@ // string representation of it. It is particularly useful since it allows for // a second parameter indicating the number of digits to which we want the // result padded with leading zeros. We will use this to write output files -// that have the form solution-XXX.gnuplot where XXX +// that have the form solution-XXX.vtu where XXX // denotes the number of the time step and always consists of three digits // even if we are still in the single or double digit time steps. #include @@ -163,90 +157,67 @@ namespace Step23 class InitialValuesU : public Function { public: - InitialValuesU() - : Function() - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point & /*p*/, + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + return 0; + } }; + template class InitialValuesV : public Function { public: - InitialValuesV() - : Function() - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point & /*p*/, + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + return 0; + } }; - template - double InitialValuesU::value(const Point & /*p*/, - const unsigned int component) const - { - (void)component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - return 0; - } - - - - template - double InitialValuesV::value(const Point & /*p*/, - const unsigned int component) const - { - (void)component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - return 0; - } - - - // Secondly, we have the right hand side forcing term. Boring as we are, we // choose zero here as well: template class RightHandSide : public Function { public: - RightHandSide() - : Function() - {} - - virtual double value(const Point & p, - const unsigned int component = 0) const override; + virtual double value(const Point & /*p*/, + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + return 0; + } }; - template - double RightHandSide::value(const Point & /*p*/, - const unsigned int component) const - { - (void)component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - return 0; - } - - - // Finally, we have boundary values for $u$ and $v$. They are as described // in the introduction, one being the time derivative of the other: template class BoundaryValuesU : public Function { public: - BoundaryValuesU() - : Function() - {} - virtual double value(const Point & p, - const unsigned int component = 0) const override; + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + + if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) && + (p[1] > -1. / 3)) + return std::sin(this->get_time() * 4 * numbers::PI); + else + return 0; + } }; @@ -255,48 +226,22 @@ namespace Step23 class BoundaryValuesV : public Function { public: - BoundaryValuesV() - : Function() - {} - virtual double value(const Point & p, - const unsigned int component = 0) const override; + const unsigned int component = 0) const override + { + (void)component; + Assert(component == 0, ExcIndexRange(component, 0, 1)); + + if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) && + (p[1] > -1. / 3)) + return (std::cos(this->get_time() * 4 * numbers::PI) * 4 * numbers::PI); + else + return 0; + } }; - template - double BoundaryValuesU::value(const Point & p, - const unsigned int component) const - { - (void)component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - - if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) && - (p[1] > -1. / 3)) - return std::sin(this->get_time() * 4 * numbers::PI); - else - return 0; - } - - - - template - double BoundaryValuesV::value(const Point & p, - const unsigned int component) const - { - (void)component; - Assert(component == 0, ExcIndexRange(component, 0, 1)); - - if ((this->get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) && - (p[1] > -1. / 3)) - return (std::cos(this->get_time() * 4 * numbers::PI) * 4 * numbers::PI); - else - return 0; - } - - - // @sect3{Implementation of the WaveEquation class} // The implementation of the actual logic is actually fairly short, since we @@ -365,9 +310,12 @@ namespace Step23 // integration. Note that in many respects these functions are better than // what we would usually do in application programs, for example because // they automatically parallelize building the matrices if multiple - // processors are available in a machine. The matrices for solving linear - // systems will be filled in the run() method because we need to re-apply - // boundary conditions every time step. + // processors are available in a machine: for more information see the + // documentation of WorkStream or the + // @ref threads "Parallel computing with multiple processors" + // module. The matrices for solving linear systems will be filled in the + // run() method because we need to re-apply boundary conditions every time + // step. mass_matrix.reinit(sparsity_pattern); laplace_matrix.reinit(sparsity_pattern); matrix_u.reinit(sparsity_pattern); @@ -396,6 +344,7 @@ namespace Step23 } + // @sect4{WaveEquation::solve_u and WaveEquation::solve_v} // The next two functions deal with solving the linear systems associated @@ -423,6 +372,7 @@ namespace Step23 } + template void WaveEquation::solve_v() { @@ -456,9 +406,17 @@ namespace Step23 data_out.build_patches(); const std::string filename = - "solution-" + Utilities::int_to_string(timestep_number, 3) + ".gnuplot"; + "solution-" + Utilities::int_to_string(timestep_number, 3) + ".vtu"; + // Like step-15, since we write output at every time step (and the system + // we have to solve is relatively easy), we instruct DataOut to use the + // zlib compression algorithm that is optimized for speed instead of disk + // usage since otherwise plotting the output becomes a bottleneck: + DataOutBase::VtkFlags vtk_flags; + vtk_flags.compression_level = + DataOutBase::VtkFlags::ZlibCompressionLevel::best_speed; + data_out.set_flags(vtk_flags); std::ofstream output(filename); - data_out.write_gnuplot(output); + data_out.write_vtu(output); } @@ -530,7 +488,7 @@ namespace Step23 RightHandSide rhs_function; rhs_function.set_time(time); VectorTools::create_right_hand_side(dof_handler, - QGauss(2), + QGauss(fe.degree + 1), rhs_function, tmp); forcing_terms = tmp; @@ -538,7 +496,7 @@ namespace Step23 rhs_function.set_time(time - time_step); VectorTools::create_right_hand_side(dof_handler, - QGauss(2), + QGauss(fe.degree + 1), rhs_function, tmp); -- 2.39.5