step-8.h

Go to the documentation of this file.

 false);
 *       sparsity_pattern.copy_from(dsp);
 *   
 *       system_matrix.reinit(sparsity_pattern);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-ElasticProblemassemble_system"></a> 
 * <h4>ElasticProblem::assemble_system</h4>
 * 
 
 * 
 * The big changes in this program are in the creation of matrix and right
 * hand side, since they are problem-dependent. We will go through that
 * process step-by-step, since it is a bit more complicated than in previous
 * examples.
 *   
 
 * 
 * The first parts of this function are the same as before, however: setting
 * up a suitable quadrature formula, initializing an FEValues object for the
 * (vector-valued) finite element we use as well as the quadrature object,
 * and declaring a number of auxiliary arrays. In addition, we declare the
 * ever same two abbreviations: <code>n_q_points</code> and
 * <code>dofs_per_cell</code>. The number of degrees of freedom per cell we
 * now obviously ask from the composed finite element rather than from the
 * underlying scalar Q1 element. Here, it is <code>dim</code> times the
 * number of degrees of freedom per cell of the Q1 element, though this is
 * not explicit knowledge we need to care about:
 * 
 * @code
 *     template <int dim>
 *     void ElasticProblem<dim>::assemble_system()
 *     {
 *       const QGauss<dim> quadrature_formula(fe.degree + 1);
 *   
 *       FEValues<dim> fe_values(fe,
 *                               quadrature_formula,
 *                               update_values | update_gradients |
 *                                 update_quadrature_points | update_JxW_values);
 *   
 *       const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
 *       const unsigned int n_q_points    = quadrature_formula.size();
 *   
 *       FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
 *       Vector<double>     cell_rhs(dofs_per_cell);
 *   
 *       std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
 *   
 * @endcode
 * 
 * As was shown in previous examples as well, we need a place where to
 * store the values of the coefficients at all the quadrature points on a
 * cell. In the present situation, we have two coefficients, lambda and
 * mu.
 * 
 * @code
 *       std::vector<double> lambda_values(n_q_points);
 *       std::vector<double> mu_values(n_q_points);
 *   
 * @endcode
 * 
 * Well, we could as well have omitted the above two arrays since we will
 * use constant coefficients for both lambda and mu, which can be declared
 * like this. They both represent functions always returning the constant
 * value 1.0. Although we could omit the respective factors in the
 * assemblage of the matrix, we use them here for purpose of
 * demonstration.
 * 
 * @code
 *       Functions::ConstantFunction<dim> lambda(1.), mu(1.);
 *   
 * @endcode
 * 
 * Like the two constant functions above, we will call the function
 * right_hand_side just once per cell to make things simpler.
 * 
 * @code
 *       std::vector<Tensor<1, dim>> rhs_values(n_q_points);
 *   
 * @endcode
 * 
 * Now we can begin with the loop over all cells:
 * 
 * @code
 *       for (const auto &cell : dof_handler.active_cell_iterators())
 *         {
 *           fe_values.reinit(cell);
 *   
 *           cell_matrix = 0;
 *           cell_rhs    = 0;
 *   
 * @endcode
 * 
 * Next we get the values of the coefficients at the quadrature
 * points. Likewise for the right hand side:
 * 
 * @code
 *           lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
 *           mu.value_list(fe_values.get_quadrature_points(), mu_values);
 *           right_hand_side(fe_values.get_quadrature_points(), rhs_values);
 *   
 * @endcode
 * 
 * Then assemble the entries of the local @ref GlossStiffnessMatrix "stiffness matrix" and right
 * hand side vector. This follows almost one-to-one the pattern
 * described in the introduction of this example.  One of the few
 * comments in place is that we can compute the number
 * <code>comp(i)</code>, i.e. the index of the only nonzero vector
 * component of shape function <code>i</code> using the
 * <code>fe.system_to_component_index(i).first</code> function call
 * below.
 *         
 
 * 
 * (By accessing the <code>first</code> variable of the return value
 * of the <code>system_to_component_index</code> function, you might
 * already have guessed that there is more in it. In fact, the
 * function returns a <code>std::pair@<unsigned int, unsigned
 * int@></code>, of which the first element is <code>comp(i)</code>
 * and the second is the value <code>base(i)</code> also noted in the
 * introduction, i.e.  the index of this shape function within all the
 * shape functions that are nonzero in this component,
 * i.e. <code>base(i)</code> in the diction of the introduction. This
 * is not a number that we are usually interested in, however.)
 *         
 
 * 
 * With this knowledge, we can assemble the local matrix
 * contributions:
 * 
 * @code
 *           for (const unsigned int i : fe_values.dof_indices())
 *             {
 *               const unsigned int component_i =
 *                 fe.system_to_component_index(i).first;
 *   
 *               for (const unsigned int j : fe_values.dof_indices())
 *                 {
 *                   const unsigned int component_j =
 *                     fe.system_to_component_index(j).first;
 *   
 *                   for (const unsigned int q_point :
 *                        fe_values.quadrature_point_indices())
 *                     {
 *                       cell_matrix(i, j) +=
 * @endcode
 * 
 * The first term is @f$(\lambda \partial_i u_i, \partial_j
 * v_j) + (\mu \partial_i u_j, \partial_j v_i)@f$. Note
 * that <code>shape_grad(i,q_point)</code> returns the
 * gradient of the only nonzero component of the i-th
 * shape function at quadrature point q_point. The
 * component <code>comp(i)</code> of the gradient, which
 * is the derivative of this only nonzero vector
 * component of the i-th shape function with respect to
 * the comp(i)th coordinate is accessed by the appended
 * brackets.
 * 
 * @code
 *                         (                                                  
 *                           (fe_values.shape_grad(i, q_point)[component_i] * 
 *                            fe_values.shape_grad(j, q_point)[component_j] * 
 *                            lambda_values[q_point])                         
 *                           +                                                
 *                           (fe_values.shape_grad(i, q_point)[component_j] * 
 *                            fe_values.shape_grad(j, q_point)[component_i] * 
 *                            mu_values[q_point])                             
 *                           +                                                
 * @endcode
 * 
 * The second term is @f$(\mu \nabla u_i, \nabla
 * v_j)@f$. We need not access a specific component of
 * the gradient, since we only have to compute the
 * scalar product of the two gradients, of which an
 * overloaded version of <tt>operator*</tt> takes
 * care, as in previous examples.
 *                         
 
 * 
 * Note that by using the <tt>?:</tt> operator, we only
 * do this if <tt>component_i</tt> equals
 * <tt>component_j</tt>, otherwise a zero is added
 * (which will be optimized away by the compiler).
 * 
 * @code
 *                           ((component_i == component_j) ?        
 *                              (fe_values.shape_grad(i, q_point) * 
 *                               fe_values.shape_grad(j, q_point) * 
 *                               mu_values[q_point]) :              
 *                              0)                                  
 *                           ) *                                    
 *                         fe_values.JxW(q_point);                  
 *                     }
 *                 }
 *             }
 *   
 * @endcode
 * 
 * Assembling the right hand side is also just as discussed in the
 * introduction:
 * 
 * @code
 *           for (const unsigned int i : fe_values.dof_indices())
 *             {
 *               const unsigned int component_i =
 *                 fe.system_to_component_index(i).first;
 *   
 *               for (const unsigned int q_point :
 *                    fe_values.quadrature_point_indices())
 *                 cell_rhs(i) += fe_values.shape_value(i, q_point) *
 *                                rhs_values[q_point][component_i] *
 *                                fe_values.JxW(q_point);
 *             }
 *   
 * @endcode
 * 
 * The transfer from local degrees of freedom into the global matrix
 * and right hand side vector does not depend on the equation under
 * consideration, and is thus the same as in all previous
 * examples.
 * 
 * @code
 *           cell->get_dof_indices(local_dof_indices);
 *           constraints.distribute_local_to_global(
 *             cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
 *         }
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-ElasticProblemsolve"></a> 
 * <h4>ElasticProblem::solve</h4>
 * 
 
 * 
 * The solver does not care about where the system of equations comes from, as
 * long as it is positive definite and symmetric (which are the
 * requirements for the use of the CG solver), which the system indeed
 * is. Therefore, we need not change anything.
 * 
 * @code
 *     template <int dim>
 *     void ElasticProblem<dim>::solve()
 *     {
 *       SolverControl            solver_control(1000, 1e-6 * system_rhs.l2_norm());
 *       SolverCG<Vector<double>> cg(solver_control);
 *   
 *       PreconditionSSOR<SparseMatrix<double>> preconditioner;
 *       preconditioner.initialize(system_matrix, 1.2);
 *   
 *       cg.solve(system_matrix, solution, system_rhs, preconditioner);
 *   
 *       constraints.distribute(solution);
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-ElasticProblemrefine_grid"></a> 
 * <h4>ElasticProblem::refine_grid</h4>
 * 
 
 * 
 * The function that does the refinement of the grid is the same as in the
 * @ref step_6 "step-6" example. The quadrature formula is adapted to the linear elements
 * again. Note that the error estimator by default adds up the estimated
 * obtained from all components of the finite element solution, i.e., it
 * uses the displacement in all directions with the same weight. If we would
 * like the grid to be adapted to the x-displacement only, we could pass the
 * function an additional parameter which tells it to do so and do not
 * consider the displacements in all other directions for the error
 * indicators. However, for the current problem, it seems appropriate to
 * consider all displacement components with equal weight.
 * 
 * @code
 *     template <int dim>
 *     void ElasticProblem<dim>::refine_grid()
 *     {
 *       Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
 *   
 *       KellyErrorEstimator<dim>::estimate(dof_handler,
 *                                          QGauss<dim - 1>(fe.degree + 1),
 *                                          {},
 *                                          solution,
 *                                          estimated_error_per_cell);
 *   
 *       GridRefinement::refine_and_coarsen_fixed_number(triangulation,
 *                                                       estimated_error_per_cell,
 *                                                       0.3,
 *                                                       0.03);
 *   
 *       triangulation.execute_coarsening_and_refinement();
 *     }
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-ElasticProblemoutput_results"></a> 
 * <h4>ElasticProblem::output_results</h4>
 * 
 
 * 
 * The output happens mostly as has been shown in previous examples
 * already. The only difference is that the solution function is vector
 * valued. The DataOut class takes care of this automatically, but we have
 * to give each component of the solution vector a different name.
 *   
 
 * 
 * To do this, the DataOut::add_vector() function wants a vector of
 * strings. Since the number of components is the same as the number
 * of dimensions we are working in, we use the <code>switch</code>
 * statement below.
 *   
 
 * 
 * We note that some graphics programs have restriction on what
 * characters are allowed in the names of variables. deal.II therefore
 * supports only the minimal subset of these characters that is supported
 * by all programs. Basically, these are letters, numbers, underscores,
 * and some other characters, but in particular no whitespace and
 * minus/hyphen. The library will throw an exception otherwise, at least
 * if in debug mode.
 *   
 
 * 
 * After listing the 1d, 2d, and 3d case, it is good style to let the
 * program die if we run into a case which we did not consider. You have
 * previously already seen the use of the `Assert` macro that generates
 * aborts the program with an error message if a condition is not satisfied
 * (see @ref step_5 "step-5", for example). We could use this in the `default` case
 * below, in the form `Assert(false, ExcNotImplemented())` -- in other words,
 * the "condition" here is always `false`, and so the assertion always fails
 * and always aborts the program whenever it gets to the default statement.
 * This is perhaps more difficult to read than necessary, and consequently
 * there is a short-cut: `DEAL_II_NOT_IMPLEMENTED()`. It does the same
 * as the form above (with the minor difference that it also aborts the
 * program in release mode). It is written in all-caps because that makes
 * it stand out visually (and also because it is not actually a function,
 * but a macro).
 * 
 * @code
 *     template <int dim>
 *     void ElasticProblem<dim>::output_results(const unsigned int cycle) const
 *     {
 *       DataOut<dim> data_out;
 *       data_out.attach_dof_handler(dof_handler);
 *   
 *       std::vector<std::string> solution_names;
 *       switch (dim)
 *         {
 *           case 1:
 *             solution_names.emplace_back("displacement");
 *             break;
 *           case 2:
 *             solution_names.emplace_back("x_displacement");
 *             solution_names.emplace_back("y_displacement");
 *             break;
 *           case 3:
 *             solution_names.emplace_back("x_displacement");
 *             solution_names.emplace_back("y_displacement");
 *             solution_names.emplace_back("z_displacement");
 *             break;
 *           default:
 *             DEAL_II_NOT_IMPLEMENTED();
 *         }
 *   
 * @endcode
 * 
 * After setting up the names for the different components of the
 * solution vector, we can add the solution vector to the list of
 * data vectors scheduled for output. Note that the following
 * function takes a vector of strings as second argument, whereas
 * the one which we have used in all previous examples accepted a
 * string there. (In fact, the function we had used before would
 * convert the single string into a vector with only one element
 * and forwards that to the other function.)
 * 
 * @code
 *       data_out.add_data_vector(solution, solution_names);
 *       data_out.build_patches();
 *   
 *       std::ofstream output("solution-" + std::to_string(cycle) + ".vtk");
 *       data_out.write_vtk(output);
 *     }
 *   
 *   
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-ElasticProblemrun"></a> 
 * <h4>ElasticProblem::run</h4>
 * 
 
 * 
 * The <code>run</code> function does the same things as in @ref step_6 "step-6", for
 * example. This time, we use the square [-1,1]^d as domain, and we refine
 * it globally four times before starting the first iteration.
 *   
 
 * 
 * The reason for refining is a bit accidental: we use the QGauss
 * quadrature formula with two points in each direction for integration of the
 * right hand side; that means that there are four quadrature points on each
 * cell (in 2d). If we only refine the initial grid once globally, then there
 * will be only four quadrature points in each direction on the
 * domain. However, the right hand side function was chosen to be rather
 * localized and in that case, by pure chance, it happens that all quadrature
 * points lie at points where the right hand side function is zero (in
 * mathematical terms, the quadrature points happen to be at points outside
 * the <i>support</i> of the right hand side function). The right hand side
 * vector computed with quadrature will then contain only zeroes (even though
 * it would of course be nonzero if we had computed the right hand side vector
 * exactly using the integral) and the solution of the system of
 * equations is the zero vector, i.e., a finite element function that is zero
 * everywhere. In a sense, we
 * should not be surprised that this is happening since we have chosen
 * an initial grid that is totally unsuitable for the problem at hand.
 *   
 
 * 
 * The unfortunate thing is that if the discrete solution is constant, then
 * the error indicators computed by the KellyErrorEstimator class are zero
 * for each cell as well, and the call to
 * Triangulation::refine_and_coarsen_fixed_number() will not flag any cells
 * for refinement (why should it if the indicated error is zero for each
 * cell?). The grid in the next iteration will therefore consist of four
 * cells only as well, and the same problem occurs again.
 *   
 
 * 
 * The conclusion needs to be: while of course we will not choose the
 * initial grid to be well-suited for the accurate solution of the problem,
 * we must at least choose it such that it has the chance to capture the
 * important features of the solution. In this case, it needs to be able to
 * see the right hand side. Thus, we refine globally four times. (Any larger
 * number of global refinement steps would of course also work.)
 * 
 * @code
 *     template <int dim>
 *     void ElasticProblem<dim>::run()
 *     {
 *       for (unsigned int cycle = 0; cycle < 8; ++cycle)
 *         {
 *           std::cout << "Cycle " << cycle << ':' << std::endl;
 *   
 *           if (cycle == 0)
 *             {
 *               GridGenerator::hyper_cube(triangulation, -1, 1);
 *               triangulation.refine_global(4);
 *             }
 *           else
 *             refine_grid();
 *   
 *           std::cout << "   Number of active cells:       "
 *                     << triangulation.n_active_cells() << std::endl;
 *   
 *           setup_system();
 *   
 *           std::cout << "   Number of degrees of freedom: " << dof_handler.n_dofs()
 *                     << std::endl;
 *   
 *           assemble_system();
 *           solve();
 *           output_results(cycle);
 *         }
 *     }
 *   } // namespace Step8
 *   
 * @endcode
 * 
 * 
 * <a name="step_8-Thecodemaincodefunction"></a> 
 * <h3>The <code>main</code> function</h3>
 * 
 
 * 
 * After closing the <code>Step8</code> namespace in the last line above, the
 * following is the main function of the program and is again exactly like in
 * @ref step_6 "step-6" (apart from the changed class names, of course).
 * 
 * @code
 *   int main()
 *   {
 *     try
 *       {
 *         Step8::ElasticProblem<2> elastic_problem_2d;
 *         elastic_problem_2d.run();
 *       }
 *     catch (std::exception &exc)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Exception on processing: " << std::endl
 *                   << exc.what() << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *   
 *         return 1;
 *       }
 *     catch (...)
 *       {
 *         std::cerr << std::endl
 *                   << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         std::cerr << "Unknown exception!" << std::endl
 *                   << "Aborting!" << std::endl
 *                   << "----------------------------------------------------"
 *                   << std::endl;
 *         return 1;
 *       }
 *   
 *     return 0;
 *   }
 * @endcode
<a name="step_8-Results"></a><h1>Results</h1>
 
 
 
There is not much to be said about the results of this program, other than
that they look nice. All images were made using VisIt from the
output files that the program wrote to disk. The first two pictures show
the @f$x@f$- and @f$y@f$-displacements as scalar components:
 
<table width="100%">
<tr>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-8.x.png" alt="">
</td>
<td>
<img src="https://www.dealii.org/images/steps/developer/step-8.y.png" alt="">
</td>
</tr>
</table>
 
 
You can clearly see the sources of @f$x@f$-displacement around @f$x=0.5@f$ and
@f$x=-0.5@f$, and of @f$y@f$-displacement at the origin.
 
What one frequently would like to do is to show the displacement as a vector
field, i.e., vectors that for each point illustrate the direction and magnitude
of displacement. Unfortunately, that's a bit more involved. To understand why
this is so, remember that we have just defined our finite element as a
collection of two  components (in <code>dim=2</code> dimensions). Nowhere have
we said that this is not just a pressure and a concentration (two scalar
quantities) but that the two components actually are the parts of a
vector-valued quantity, namely the displacement. Absent this knowledge, the
DataOut class assumes that all individual variables we print are separate
scalars, and VisIt and Paraview then faithfully assume that this is indeed what it is. In
other words, once we have written the data as scalars, there is nothing in
these programs that allows us to paste these two scalar fields back together as a
vector field. Where we would have to attack this problem is at the root,
namely in <code>ElasticProblem::output_results()</code>. We won't do so here but
instead refer the reader to the @ref step_22 "step-22" program where we show how to do this
for a more general situation. That said, we couldn't help generating the data
anyway that would show how this would look if implemented as discussed in
@ref step_22 "step-22". The vector field then looks like this (VisIt and Paraview
randomly select a few
hundred vertices from which to draw the vectors; drawing them from each
individual vertex would make the picture unreadable):
 
<img src="https://www.dealii.org/images/steps/developer/step-8.vectors.png" alt="">
 
 
We note that one may have intuitively expected the
solution to be symmetric about the @f$x@f$- and @f$y@f$-axes since the @f$x@f$- and
@f$y@f$-forces are symmetric with respect to these axes. However, the force
considered as a vector is not symmetric and consequently neither is
the solution.
 *
 *
<a name="step_8-PlainProg"></a>
<h1> The plain program</h1>
@include "step-8.cc"
*/