simulating problems in which the deformation becomes <i>large</i>.
@note The model we will consider below is not founded on anything that
-would be mathematical sound: we will consider a model in which we
+would be mathematically sound: we will consider a model in which we
produce a small deformation, deform the physical coordinates of the
body by this deformation, and then consider the next loading step
again as a linear problem. This isn't consistent, since the assumption
Converter from deal.II intermediate format to other graphics formats.
-Usage: ./\step-19 [-p parameter_file] list_of_input_files [-x output_format] output_file
+Usage: ./\step-19 [-p parameter_file] list_of_input_files [-x output_format] -o output_file
-examples/\step-18> ../\step-19/\step-19 solution-0001.0000.d2 -x gmv solution-0001.0000.gmv
-examples/\step-18> ../\step-19/\step-19 solution-0002.0000.d2 -x gmv solution-0002.0000.gmv
+examples/\step-18> ../\step-19/\step-19 solution-0001.0000.d2 -x gmv -o solution-0001.0000.gmv
+examples/\step-18> ../\step-19/\step-19 solution-0002.0000.d2 -x gmv -o solution-0002.0000.gmv
[...]
@endcode
Of course, since we have run the program only in sequential mode, we
// in the form $C_{ijkl} = \mu (\delta_{ik} \delta_{jl} + \delta_{il}
// \delta_{jk}) + \lambda \delta_{ij} \delta_{kl}$. This tensor maps
// symmetric tensor of rank 2 to symmetric tensors of rank 2. A function
- // implementing its creation for given values of the Lame constants lambda
- // and mu is straightforward:
+ // implementing its creation for given values of the Lame constants $\lambda$
+ // and $\mu$ is straightforward:
template <int dim>
SymmetricTensor<4,dim>
get_stress_strain_tensor (const double lambda, const double mu)
// One difference of this program is that we declare the quadrature
// formula in the class declaration. The reason is that in all the other
// programs, it didn't do much harm if we had used different quadrature
- // formulas when computing the matrix and the righ hand side, for
+ // formulas when computing the matrix and the right hand side, for
// example. However, in the present case it does: we store information in
// the quadrature points, so we have to make sure all parts of the program
// agree on where they are and how many there are on each cell. Thus, let
// As a last piece of data, let us also add the partitioning of the domain
// into subdomains associated with the processors if this is a parallel
// job. This works in the exact same way as in the step-17 program:
- std::vector<unsigned int> partition_int (triangulation.n_active_cells());
+ std::vector<types::subdomain_id> partition_int (triangulation.n_active_cells());
GridTools::get_subdomain_association (triangulation, partition_int);
const Vector<double> partitioning(partition_int.begin(),
partition_int.end());
prm.get_bool ("Dummy generate output");
}
prm.leave_subsection ();
- // We would assign the result to a variable, or course, but don't
+ // We would assign the result to a variable, of course, but don't
// here in order not to generate an unused variable that the
// compiler might warn about.
//
down time stepping schemes. Note that we do not have boundary
conditions for $v$ at first. However, we could enforce $v=\frac{\partial
g}{\partial t}$ on the boundary. It turns out in numerical examples that this
-as actually necessary: without doing so the solution doesn't look particularly
+is actually necessary: without doing so the solution doesn't look particularly
wrong, but the Crank-Nicolson scheme does not conserve energy if one doesn't
enforce these boundary conditions.
\\
A^{n,n-1}_{ij} &=& (\nabla\phi_i^n, \nabla\phi_j^{n-1}),
\\
- F^n_{ij} &=& (f^n,\phi_i^n),
+ F^n_{i} &=& (f^n,\phi_i^n),
\\
- F^{n-1}_{ij} &=& (f^{n-1},\phi_i^n).
+ F^{n-1}_{i} &=& (f^{n-1},\phi_i^n).
@f}
If we solve these two equations, we can move the solution one step
To make such a change, we would have to compute the mass and Laplace
matrices with a variable coefficient. Fortunately, this isn't too hard: the
- functions MatrixTools::create_laplace_matrix and
- MatrixTools::create_vector_tools have additional default parameters that can
+ functions MatrixCreator::create_laplace_matrix and
+ MatrixCreator::create_mass_matrix have additional default parameters that can
be used to pass non-constant coefficient functions to them. The required
changes are therefore relatively small. On the other hand, care must be
taken again to make sure the time step is within the allowed range.
// Here are the only three include files of some new interest: The first one
// is already used, for example, for the
// VectorTools::interpolate_boundary_values and
-// VectorTools::apply_boundary_values functions. However, we here use another
+// MatrixTools::apply_boundary_values functions. However, we here use another
// function in that class, VectorTools::project to compute our initial values
// as the $L^2$ projection of the continuous initial values. Furthermore, we
// use VectorTools::create_right_hand_side to generate the integrals
// mass and Laplace matrices, although it would have only taken copying the
// relevant code from any number of previous tutorial programs. Rather, we
// want to focus on the things that are truly new to this program and
-// therefore use the MatrixTools::create_mass_matrix and
-// MatrixTools::create_laplace_matrix functions. They are declared here:
+// therefore use the MatrixCreator::create_mass_matrix and
+// MatrixCreator::create_laplace_matrix functions. They are declared here:
#include <deal.II/numerics/matrix_tools.h>
// Finally, here is an include file that contains all sorts of tool functions
v = \frac{\partial\bar{p}}{\partial t}
@f]
-With the second variables, one then transform the forward problem into
+With the second variable, one then transforms the forward problem into
two separate equations:
@f{eqnarray*}
\bar{p}_{t} - v & = & 0 \\
// pressure. In the physical setting considered in this program, these are
// small absorbers, which we model as a series of little circles where we
// assume that the pressure surplus is one, whereas no absorption and
- // therefore no pressure surplus is anywhere else. This is how we do things
+ // therefore no pressure surplus is everywhere else. This is how we do things
// (note that if we wanted to expand this program to not only compile but
// also to run, we would have to initialize the sources with
// three-dimensional source locations):
//
// A first observation would be that this matrix is much sparser than the
// regular mass matrix, since none of the shape functions with purely
- // interior support contributes to this matrix. We could therefore
+ // interior support contribute to this matrix. We could therefore
// optimize the storage pattern to this situation and build up a second
// sparsity pattern that only contains the nonzero entries that we
// need. There is a trade-off to make here: first, we would have to have a
// consumption of the program, the introduction of a few temporary vectors
// isn't doing much harm.
//
- // The only changes to this function are: First, that we do not have to
+ // The only changes to this function are: first, that we do not have to
// project initial values for the velocity $v$, since we know that it is
// zero. And second that we evaluate the solution at the detector locations
// computed in the constructor. This is done using the
\mbox{ Find } \delta u^n_l \mbox{ s.t. } F'(u^n_l)\delta u^n_l = -F(u^n_l)
\mbox{, set } u^n_{l+1} = u^n_l + \delta u^n_l.
\f}
-The iteration can be initialized with the old time step, i.e. $u^{n,0} =
-u^{n-1}$, and eventually it will produce a solution to the first equation of
+The iteration can be initialized with the old time step, i.e. $u^n_0 = u^{n-1}$,
+and eventually it will produce a solution to the first equation of
the split formulation (see above). For the time discretizaion of the
sine-Gordon equation under consideration here, we have that
\f{eqnarray*}
&\mbox{ Find}& \delta u^n_l \in H^1(\Omega) \mbox{ s.t. }
\left( F'(u^n_l)\delta u^n_l, \varphi \right)_{\Omega}
= -\left(F(u^n_l), \varphi \right)_{\Omega} \;\forall\varphi\in H^1(\Omega),
- \mbox{ set } u^n_{l+1} = u^n_l + \delta u^n_l,\; u^{n,0} = u^{n-1}.\\
+ \mbox{ set } u^n_{l+1} = u^n_l + \delta u^n_l,\; u^n_0 = u^{n-1}.\\
&\mbox{ Find}& v^n \in H^1(\Omega) \mbox{ s.t. }
\left( v^n, \varphi \right)_{\Omega} = \left( v^{n-1}, \varphi \right)_{\Omega}
- k\theta\left( \nabla u^n, \nabla\varphi \right)_{\Omega}
\f{eqnarray*}
F_h(U^{n,l}) &=& \left[ M+k^2\theta^2A \right] U^{n,l} -
\left[ M-k^2\theta(1-\theta)A \right] U^{n-1} - k MV^{n-1}
- + k^2\theta S(u^{n,l}, u^{n-1}),\\
+ + k^2\theta S(u^n_l, u^{n-1}),\\
F_h'(U^{n,l}) &=& M+k^2\theta^2A
- - k^2\theta^2N(u^{n,l},u^{n-1})
+ - k^2\theta^2N(u^n_l,u^{n-1})
\f}
Again, note that the first matrix equation above is, in fact, the
defition of an iterative procedure, so it is solved multiple times
What solvers can we use for the first equation? Let's look at the matrix we
have to invert:
@f[
- (M-k^2\theta^2N)_{ij} =
+ (M+k^2\theta^2(A-N))_{ij} =
\int_\Omega (1-k^2\theta^2 \cos \alpha)
- \varphi_i\varphi_j \; dx,
+ \varphi_i\varphi_j \; dx+\int_\Omega \nabla\varphi_i\nabla\varphi_j \; dx,
@f]
for some $\alpha$ that depends on the present and previous solution. First,
note that the matrix is symmetric. In addition, if the time step $k$ is small
// @sect4{SineGordonProblem::assemble_system}
- // This functions assembles the system matrix and right-hand side vector for
+ // This function assembles the system matrix and right-hand side vector for
// each iteration of Newton's method. The reader should refer to the
// Introduction for the explicit formulas for the system matrix and
// right-hand side.
{
// First we assemble the Jacobian matrix $F'_h(U^{n,l})$, where $U^{n,l}$
// is stored in the vector <code>solution</code> for convenience.
- system_matrix = 0;
system_matrix.copy_from (mass_matrix);
system_matrix.add (std::pow(time_step*theta,2), laplace_matrix);
// Then, we compute the right-hand side vector $-F_h(U^{n,l})$.
system_rhs = 0;
- tmp_matrix = 0;
tmp_matrix.copy_from (mass_matrix);
tmp_matrix.add (std::pow(time_step*theta,2), laplace_matrix);
tmp_matrix.vmult (tmp_vector, solution);
system_rhs += tmp_vector;
- tmp_matrix = 0;
tmp_matrix.copy_from (mass_matrix);
tmp_matrix.add (-std::pow(time_step,2)*theta*(1-theta), laplace_matrix);
- tmp_vector = 0;
tmp_matrix.vmult (tmp_vector, old_solution);
system_rhs -= tmp_vector;
system_rhs.add (-time_step, M_x_velocity);
- tmp_vector = 0;
compute_nl_term (old_solution, solution, tmp_vector);
system_rhs.add (std::pow(time_step,2)*theta, tmp_vector);
const Vector<double> &new_data,
Vector<double> &nl_term) const
{
+ nl_term = 0;
const QGauss<dim> quadrature_formula (3);
FEValues<dim> fe_values (fe, quadrature_formula,
update_values |
for (; cell!=endc; ++cell)
{
+ local_nl_term = 0;
// Once we re-initialize our <code>FEValues</code> instantiation to
// the current cell, we make use of the
// <code>get_function_values</code> routine to get the values of the
for (unsigned int i=0; i<dofs_per_cell; ++i)
nl_term(local_dof_indices[i]) += local_nl_term(i);
-
- local_nl_term = 0;
}
}
for (; cell!=endc; ++cell)
{
+ local_nl_matrix = 0;
// Again, first we re-initialize our <code>FEValues</code>
// instantiation to the current cell.
fe_values.reinit (cell);
for (unsigned int j=0; j<dofs_per_cell; ++j)
nl_matrix.add(local_dof_indices[i], local_dof_indices[j],
local_nl_matrix(i,j));
-
- local_nl_matrix = 0;
}
}
PreconditionSSOR<> preconditioner;
preconditioner.initialize(system_matrix, 1.2);
- solution_update = 0;
cg.solve (system_matrix, solution_update,
system_rhs,
preconditioner);
laplace_matrix.vmult (tmp_vector, solution);
M_x_velocity.add (-time_step*theta, tmp_vector);
- tmp_vector = 0;
laplace_matrix.vmult (tmp_vector, old_solution);
M_x_velocity.add (-time_step*(1-theta), tmp_vector);
- tmp_vector = 0;
compute_nl_term (old_solution, solution, tmp_vector);
M_x_velocity.add (-time_step, tmp_vector);
[For those who wonder how this can be achieved in practice without
inadvertently getting slightly larger than one and triggering a nuclear bomb:
first, fission processes happen on different time scales. While most neutrons
-are releases very quickly after a fission event, a small number of neutrons
+are released very quickly after a fission event, a small number of neutrons
are only released by daughter nuclei after several further decays, up to 10-60
seconds after the fission was initiated. If one is therefore slightly beyond
$k_{\mathrm{eff}}=1$, one therefore has many seconds to react until all the
<li> If the cell $K$ is active on mesh $g$, but not $g'$, then the
basis functions $\varphi_{g'}^j$ are only defined either on the children
$K_c,0\le c<2^{\texttt{dim}}$, or on children of these children if cell $K$
- is refined more than once more on mesh $g'$.
+ is refined more than once on mesh $g'$.
Let us assume for a second that $K$ is only once more refined on mesh $g'$
than on mesh $g$. Using the fact that we use embedded finite element spaces
The original purpose of this program is to simulate the focussing properties
of an ultrasound wave generated by a transducer lens with variable
geometry. Recent applications in medical imaging use ultrasound waves not only
-for imaging porposes, but also to excite certain local effects in a
+for imaging purposes, but also to excite certain local effects in a
material, like changes in optical properties, that can then be measured by
other imaging techniques. A vital ingredient for these methods is the ability
to focus the intensity of the ultrasound wave in a particular part of the
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.
+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
ExcDimensionMismatch (computed_quantities[i].size(), 1));
Assert(uh[i].size() == 2, ExcDimensionMismatch (uh[i].size(), 2));
- computed_quantities[i](0) = sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
+ computed_quantities[i](0) = std::sqrt(uh[i](0)*uh[i](0) + uh[i](1)*uh[i](1));
}
}
projects / dealii.git / commitdiff