662 *
const unsigned int component = 0)
const override
665 *
Assert(component == 0, ExcIndexRange(component, 0, 1));
724 *
if ((this->
get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
744 *
if ((this->
get_time() <= 0.5) && (p[0] < 0) && (p[1] < 1. / 3) &&
757 * <a name=
"step_23-ImplementationofthecodeWaveEquationcodeclass"></a>
771 *
Let's start with the constructor (for an explanation of the choice of
772 * time step, see the section on Courant, Friedrichs, and Lewy in the
777 * WaveEquation<dim>::WaveEquation()
779 * , dof_handler(triangulation)
780 * , time_step(1. / 64)
782 * , timestep_number(1)
790 * <a name="step_23-WaveEquationsetup_system"></a>
791 * <h4>WaveEquation::setup_system</h4>
795 * The next function is the one that sets up the mesh, DoFHandler, and
796 * matrices and vectors at the beginning of the program, i.e. before the
797 * first time step. The first few lines are pretty much standard if you've
807 *
std::cout <<
"Number of active cells: " <<
triangulation.n_active_cells()
810 *
dof_handler.distribute_dofs(fe);
812 *
std::cout <<
"Number of degrees of freedom: " << dof_handler.n_dofs()
818 *
sparsity_pattern.copy_from(
dsp);
848 * @
ref threads
"Parallel computing with multiple processors"
854 *
mass_matrix.reinit(sparsity_pattern);
856 *
matrix_u.reinit(sparsity_pattern);
857 *
matrix_v.reinit(sparsity_pattern);
860 *
QGauss<dim>(fe.degree + 1),
863 *
QGauss<dim>(fe.degree + 1),
882 *
constraints.close();
919 *
std::cout <<
" u-equation: " << solver_control.last_step()
920 *
<<
" CG iterations." << std::endl;
933 *
std::cout <<
" v-equation: " << solver_control.last_step()
934 *
<<
" CG iterations." << std::endl;
942 * <a name=
"step_23-WaveEquationoutput_results"></a>
948 * before. The only thing worth mentioning is how here we generate a string
949 * representation of the time step number padded with leading zeros to 3
950 * character length using the Utilities::int_to_string function's
second
959 *
data_out.attach_dof_handler(dof_handler);
963 *
data_out.build_patches();
977 *
data_out.set_flags(vtk_flags);
979 *
data_out.write_vtu(output);
987 * <a name=
"step_23-WaveEquationrun"></a>
988 * <
h4>WaveEquation::run</
h4>
999 * be any simpler than that:
1002 * template <int dim>
1003 * void WaveEquation<dim>::run()
1007 * VectorTools::project(dof_handler,
1009 * QGauss<dim>(fe.degree + 1),
1010 * InitialValuesU<dim>(),
1012 * VectorTools::project(dof_handler,
1014 * QGauss<dim>(fe.degree + 1),
1015 * InitialValuesV<dim>(),
1020 * The next thing is to loop over all the time steps until we reach the
1021 * end time (@f$T=5@f$ in this case). In each time step, we first have to
1022 * solve for @f$U^n@f$, using the equation @f$(M^n + k^2\theta^2 A^n)U^n =@f$
1023 * @f$(M^{n,n-1} - k^2\theta(1-\theta) A^{n,n-1})U^{n-1} + kM^{n,n-1}V^{n-1}
1024 * +@f$ @f$k\theta \left[k \theta F^n + k(1-\theta) F^{n-1} \right]@f$. Note
1025 * that we use the same mesh for all time steps, so that @f$M^n=M^{n,n-1}=M@f$
1026 * and @f$A^n=A^{n,n-1}=A@f$. What we therefore have to do first is to add up
1027 * @f$MU^{n-1} - k^2\theta(1-\theta) AU^{n-1} + kMV^{n-1}@f$ and the forcing
1028 * terms, and put the result into the <code>system_rhs</code> vector. (For
1029 * these additions, we need a temporary vector that we declare before the
1030 * loop to avoid repeated memory allocations in each time step.)
1034 * The one thing to realize here is how we communicate the time variable
1035 * to the object describing the right hand side: each object derived from
1036 * the Function class has a time field that can be set using the
1037 * Function::set_time and read by Function::get_time. In essence, using
1038 * this mechanism, all functions of space and time are therefore
1039 * considered functions of space evaluated at a particular time. This
1040 * matches well what we typically need in finite element programs, where
1041 * we almost always work on a single time step at a time, and where it
1042 * never happens that, for example, one would like to evaluate a
1043 * space-time function for all times at any given spatial location.
1046 * Vector<double> tmp(solution_u.size());
1047 * Vector<double> forcing_terms(solution_u.size());
1049 * for (; time <= 5; time += time_step, ++timestep_number)
1051 * std::cout << "Time step " << timestep_number << " at t=" << time
1054 * mass_matrix.vmult(system_rhs, old_solution_u);
1056 * mass_matrix.vmult(tmp, old_solution_v);
1057 * system_rhs.add(time_step, tmp);
1059 * laplace_matrix.vmult(tmp, old_solution_u);
1060 * system_rhs.add(-theta * (1 - theta) * time_step * time_step, tmp);
1062 * RightHandSide<dim> rhs_function;
1063 * rhs_function.set_time(time);
1064 * VectorTools::create_right_hand_side(dof_handler,
1065 * QGauss<dim>(fe.degree + 1),
1068 * forcing_terms = tmp;
1069 * forcing_terms *= theta * time_step;
1071 * rhs_function.set_time(time - time_step);
1072 * VectorTools::create_right_hand_side(dof_handler,
1073 * QGauss<dim>(fe.degree + 1),
1077 * forcing_terms.add((1 - theta) * time_step, tmp);
1079 * system_rhs.add(theta * time_step, forcing_terms);
1083 * After so constructing the right hand side vector of the first
1084 * equation, all we have to do is apply the correct boundary
1085 * values. As for the right hand side, this is a space-time function
1086 * evaluated at a particular time, which we interpolate at boundary
1087 * nodes and then use the result to apply boundary values as we
1088 * usually do. The result is then handed off to the solve_u()
1093 * BoundaryValuesU<dim> boundary_values_u_function;
1094 * boundary_values_u_function.set_time(time);
1096 * std::map<types::global_dof_index, double> boundary_values;
1097 * VectorTools::interpolate_boundary_values(dof_handler,
1099 * boundary_values_u_function,
1104 * The matrix for solve_u() is the same in every time steps, so one
1105 * could think that it is enough to do this only once at the
1106 * beginning of the simulation. However, since we need to apply
1107 * boundary values to the linear system (which eliminate some matrix
1108 * rows and columns and give contributions to the right hand side),
1109 * we have to refill the matrix in every time steps before we
1110 * actually apply boundary data. The actual content is very simple:
1111 * it is the sum of the mass matrix and a weighted Laplace matrix:
1114 * matrix_u.copy_from(mass_matrix);
1115 * matrix_u.add(theta * theta * time_step * time_step, laplace_matrix);
1116 * MatrixTools::apply_boundary_values(boundary_values,
1126 * The second step, i.e. solving for @f$V^n@f$, works similarly, except
1127 * that this time the matrix on the left is the mass matrix (which we
1128 * copy again in order to be able to apply boundary conditions, and
1129 * the right hand side is @f$MV^{n-1} - k\left[ \theta A U^n +
1130 * (1-\theta) AU^{n-1}\right]@f$ plus forcing terms. Boundary values
1131 * are applied in the same way as before, except that now we have to
1132 * use the BoundaryValuesV class:
1135 * laplace_matrix.vmult(system_rhs, solution_u);
1136 * system_rhs *= -theta * time_step;
1138 * mass_matrix.vmult(tmp, old_solution_v);
1139 * system_rhs += tmp;
1141 * laplace_matrix.vmult(tmp, old_solution_u);
1142 * system_rhs.add(-time_step * (1 - theta), tmp);
1144 * system_rhs += forcing_terms;
1147 * BoundaryValuesV<dim> boundary_values_v_function;
1148 * boundary_values_v_function.set_time(time);
1150 * std::map<types::global_dof_index, double> boundary_values;
1151 * VectorTools::interpolate_boundary_values(dof_handler,
1153 * boundary_values_v_function,
1155 * matrix_v.copy_from(mass_matrix);
1156 * MatrixTools::apply_boundary_values(boundary_values,
1165 * Finally, after both solution components have been computed, we
1166 * output the result, compute the energy in the solution, and go on to
1167 * the next time step after shifting the present solution into the
1168 * vectors that hold the solution at the previous time step. Note the
1169 * function SparseMatrix::matrix_norm_square that can compute
1170 * @f$\left<V^n,MV^n\right>@f$ and @f$\left<U^n,AU^n\right>@f$ in one step,
1171 * saving us the expense of a temporary vector and several lines of
1177 * std::cout << " Total energy: "
1178 * << (mass_matrix.matrix_norm_square(solution_v) +
1179 * laplace_matrix.matrix_norm_square(solution_u)) /
1183 * old_solution_u = solution_u;
1184 * old_solution_v = solution_v;
1187 * } // namespace Step23
1193 * <a name="step_23-Thecodemaincodefunction"></a>
1194 * <h3>The <code>main</code> function</h3>
1198 * What remains is the main function of the program. There is nothing here
1206 *
using namespace Step23;
1211 *
catch (std::exception &exc)
1213 *
std::cerr << std::endl
1215 *
<<
"----------------------------------------------------"
1217 *
std::cerr <<
"Exception on processing: " << std::endl
1218 *
<< exc.what() << std::endl
1219 *
<<
"Aborting!" << std::endl
1220 *
<<
"----------------------------------------------------"
1227 *
std::cerr << std::endl
1229 *
<<
"----------------------------------------------------"
1231 *
std::cerr <<
"Unknown exception!" << std::endl
1232 *
<<
"Aborting!" << std::endl
1233 *
<<
"----------------------------------------------------"
1249Time step 1 at t=0.015625
1252 Total energy: 1.17887
1253Time step 2 at t=0.03125
1256 Total energy: 2.9655
1257Time step 3 at t=0.046875
1260 Total energy: 4.33761
1261Time step 4 at t=0.0625
1264 Total energy: 5.35499
1265Time step 5 at t=0.078125
1268 Total energy: 6.18652
1269Time step 6 at t=0.09375
1272 Total energy: 6.6799
1276Time step 31 at t=0.484375
1279 Total energy: 21.9068
1280Time step 32 at t=0.5
1283 Total energy: 23.3394
1284Time step 33 at t=0.515625
1287 Total energy: 23.1019
1291Time step 319 at t=4.98438
1294 Total energy: 23.1019
1298 Total energy: 23.1019
1309<
img src=
"https://www.dealii.org/images/steps/developer/step-23.movie.gif" alt=
"Animation of the solution of step 23.">
1317<a name=
"step-23-extensions"></a>
1318<a name=
"step_23-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
1332 in @
ref step_24
"step-24".
1348 functions MatrixCreator::create_laplace_matrix and
1349 MatrixCreator::create_mass_matrix have additional default parameters that can
1350 be used to pass non-constant coefficient functions to them. The required
1351 changes are therefore relatively small. On the other hand, care must be
1352 taken again to make sure the time step is within the allowed range.
1354 <li>In the in-code comments, we discussed the fact that the matrices for
1355 solving for @f$U^n@f$ and @f$V^n@f$ need to be reset in every time because of
1356 boundary conditions, even though the actual content does not change. It is
1357 possible to avoid copying by not eliminating columns in the linear systems,
1358 which is implemented by appending a @p false argument to the call:
1360 MatrixTools::apply_boundary_values(boundary_values,
1367 <li>deal.II being a library that supports adaptive meshes it would of course be
1368 nice if this program supported change the mesh every few time steps. Given the
1369 structure of the solution — a wave that travels through the domain —
1370 it would seem appropriate if we only refined the mesh where the wave currently is,
1371 and not simply everywhere. It is intuitively clear that we should be able to
1372 save a significant amount of cells this way. (Though upon further thought one
1373 realizes that this is really only the case in the initial stages of the simulation.
1374 After some time, for wave phenomena, the domain is filled with reflections of
1375 the initial wave going in every direction and filling every corner of the domain.
1376 At this point, there is in general little one can gain using local mesh
1379 To make adaptively changing meshes possible, there are basically two routes.
1380 The "correct" way would be to go back to the weak form we get using Rothe's
1435 method in @
ref step_31
"step-31",
for example.
1439<a name=
"step_23-PlainProg"></a>
#define Assert(cond, exc)
static ::ExceptionBase & ExcIndexRange(std::size_t arg1, std::size_t arg2, std::size_t arg3)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorType::value_type * end(VectorType &V)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
DEAL_II_HOST constexpr SymmetricTensor< 2, dim, Number > invert(const SymmetricTensor< 2, dim, Number > &)