769 * sparsity_pattern.copy_from(dsp);
771 * system_matrix.reinit(sparsity_pattern);
778 * <a name=
"step_8-ElasticProblemassemble_system"></a>
779 * <h4>ElasticProblem::assemble_system</h4>
783 * The big changes in
this program are in the creation of
matrix and right
784 * hand side, since they are problem-dependent. We will go through that
785 * process step-by-step, since it is a bit more complicated than in previous
790 * The
first parts of
this function are the same as before, however: setting
791 * up a suitable quadrature formula, initializing an
FEValues object for the
792 * (vector-valued) finite element we use as well as the quadrature
object,
793 * and declaring a number of auxiliary arrays. In addition, we declare the
794 * ever same two abbreviations: <code>n_q_points</code> and
795 * <code>dofs_per_cell</code>. The number of degrees of freedom per cell we
796 * now obviously ask from the composed finite element rather than from the
797 * underlying
scalar Q1 element. Here, it is <code>dim</code> times the
798 * number of degrees of freedom per cell of the Q1 element, though this is
799 * not explicit knowledge we need to care about:
803 * void ElasticProblem<dim>::assemble_system()
805 * const
QGauss<dim> quadrature_formula(fe.degree + 1);
808 * quadrature_formula,
812 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
813 *
const unsigned int n_q_points = quadrature_formula.size();
818 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
822 * As was shown in previous examples as well, we need a place where to
823 * store the
values of the coefficients at all the quadrature points on a
824 * cell. In the present situation, we have two coefficients,
lambda and
828 * std::vector<double> lambda_values(n_q_points);
829 * std::vector<double> mu_values(n_q_points);
833 * Well, we could as well have omitted the above two arrays since we will
834 * use
constant coefficients
for both
lambda and mu, which can be declared
836 *
value 1.0. Although we could omit the respective factors in the
837 * assemblage of the
matrix, we use them here
for purpose of
846 * right_hand_side just once per cell to make things simpler.
849 * std::vector<Tensor<1, dim>> rhs_values(n_q_points);
853 * Now we can
begin with the
loop over all cells:
856 *
for (
const auto &cell : dof_handler.active_cell_iterators())
865 * Next we get the
values of the coefficients at the quadrature
866 * points. Likewise
for the right hand side:
869 *
lambda.value_list(fe_values.get_quadrature_points(), lambda_values);
870 * mu.value_list(fe_values.get_quadrature_points(), mu_values);
871 * right_hand_side(fe_values.get_quadrature_points(), rhs_values);
875 * Then
assemble the entries of the local @ref GlossStiffnessMatrix
"stiffness matrix" and right
876 * hand side vector. This follows almost one-to-one the pattern
877 * described in the introduction of
this example. One of the few
878 * comments in place is that we can compute the number
879 * <code>comp(i)</code>, i.e. the
index of the only
nonzero vector
880 * component of shape function <code>i</code>
using the
881 * <code>fe.system_to_component_index(i).first</code> function call
886 * (By accessing the <code>
first</code> variable of the
return value
887 * of the <code>system_to_component_index</code> function, you might
888 * already have guessed that there is more in it. In fact, the
889 * function returns a <code>std::pair@<
unsigned int,
unsigned
890 *
int@></code>, of which the
first element is <code>comp(i)</code>
891 * and the
second is the value <code>base(i)</code> also noted in the
892 * introduction, i.e. the
index of
this shape function within all the
894 * i.e. <code>base(i)</code> in the diction of the introduction. This
895 * is not a number that we are usually interested in, however.)
903 *
for (
const unsigned int i : fe_values.dof_indices())
905 * const unsigned
int component_i =
906 * fe.system_to_component_index(i).
first;
908 *
for (
const unsigned int j : fe_values.dof_indices())
910 * const unsigned
int component_j =
911 * fe.system_to_component_index(j).
first;
913 *
for (
const unsigned int q_point :
914 * fe_values.quadrature_point_indices())
919 * The
first term is @f$(\
lambda \partial_i u_i, \partial_j
920 * v_j) + (\mu \partial_i u_j, \partial_j v_i)@f$. Note
921 * that <code>shape_grad(i,q_point)</code> returns the
923 * shape function at quadrature
point q_point. The
924 * component <code>comp(i)</code> of the
gradient, which
925 * is the derivative of this only
nonzero vector
926 * component of the i-th shape function with respect to
927 * the comp(i)th coordinate is accessed by the appended
932 * (fe_values.shape_grad(i, q_point)[component_i] *
933 * fe_values.shape_grad(j, q_point)[component_j] *
934 * lambda_values[q_point])
936 * (fe_values.shape_grad(i, q_point)[component_j] *
937 * fe_values.shape_grad(j, q_point)[component_i] *
938 * mu_values[q_point])
942 * The
second term is @f$(\mu \nabla u_i, \nabla
943 * v_j)@f$. We need not access a specific component of
944 * the
gradient, since we only have to compute the
946 * overloaded version of <tt>operator*</tt> takes
947 * care, as in previous examples.
951 * Note that by using the <tt>?:</tt> operator, we only
952 * do this if <tt>component_i</tt> equals
953 * <tt>component_j</tt>, otherwise a zero is added
954 * (which will be optimized away by the compiler).
957 * ((component_i == component_j) ?
958 * (fe_values.shape_grad(i, q_point) *
959 * fe_values.shape_grad(j, q_point) *
960 * mu_values[q_point]) :
963 * fe_values.JxW(q_point);
970 * Assembling the right hand side is also just as discussed in the
974 *
for (
const unsigned int i : fe_values.dof_indices())
976 * const unsigned
int component_i =
977 * fe.system_to_component_index(i).
first;
979 *
for (
const unsigned int q_point :
980 * fe_values.quadrature_point_indices())
981 * cell_rhs(i) += fe_values.shape_value(i, q_point) *
982 * rhs_values[q_point][component_i] *
983 * fe_values.JxW(q_point);
988 * The transfer from local degrees of freedom into the global
matrix
989 * and right hand side vector does not depend on the equation under
990 * consideration, and is thus the same as in all previous
994 * cell->get_dof_indices(local_dof_indices);
995 * constraints.distribute_local_to_global(
996 * cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1005 * <a name=
"step_8-ElasticProblemsolve"></a>
1006 * <h4>ElasticProblem::solve</h4>
1010 * The solver does not care about where the system of equations comes from, as
1012 * requirements
for the use of the CG solver), which the system indeed
1013 * is. Therefore, we need not change anything.
1016 *
template <
int dim>
1017 *
void ElasticProblem<dim>::solve()
1019 *
SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
1023 * preconditioner.
initialize(system_matrix, 1.2);
1025 * cg.solve(system_matrix, solution, system_rhs, preconditioner);
1027 * constraints.distribute(solution);
1034 * <a name=
"step_8-ElasticProblemrefine_grid"></a>
1035 * <h4>ElasticProblem::refine_grid</h4>
1039 * The function that does the refinement of the grid is the same as in the
1040 * @ref step_6
"step-6" example. The quadrature formula is adapted to the linear elements
1041 * again. Note that the error estimator by
default adds up the estimated
1042 * obtained from all components of the finite element solution, i.e., it
1043 * uses the displacement in all directions with the same weight. If we would
1044 * like the grid to be adapted to the x-displacement only, we could pass the
1045 * function an additional parameter which tells it to
do so and
do not
1046 * consider the displacements in all other directions
for the error
1047 * indicators. However,
for the current problem, it seems appropriate to
1048 * consider all displacement components with
equal weight.
1051 *
template <
int dim>
1052 *
void ElasticProblem<dim>::refine_grid()
1060 * estimated_error_per_cell);
1063 * estimated_error_per_cell,
1074 * <a name=
"step_8-ElasticProblemoutput_results"></a>
1075 * <h4>ElasticProblem::output_results</h4>
1079 * The output happens mostly as has been shown in previous examples
1080 * already. The only difference is that the solution function is vector
1081 * valued. The
DataOut class takes care of this automatically, but we have
1082 * to give each component of the solution vector a different name.
1086 * To
do this, the DataOut::add_vector() function wants a vector of
1087 * strings. Since the number of components is the same as the number
1088 * of dimensions we are working in, we use the <code>switch</code>
1093 * We note that some graphics programs have restriction on what
1094 * characters are allowed in the names of variables. deal.II therefore
1095 * supports only the minimal subset of these characters that is supported
1096 * by all programs. Basically, these are letters,
numbers, underscores,
1097 * and some other characters, but in particular no whitespace and
1098 * minus/hyphen. The library will throw an exception otherwise, at least
1103 * After listing the 1d, 2d, and 3d case, it is good style to let the
1104 * program die if we run into a case which we did not consider. You have
1105 * previously already seen the use of the `
Assert` macro that generates
1106 * aborts the program with an error message if a condition is not satisfied
1107 * (see @ref step_5 "step-5", for example). We could use this in the `default` case
1108 * below, in the form `
Assert(false, ExcNotImplemented())` -- in other words,
1109 * the "condition" here is always `false`, and so the assertion always fails
1110 * and always aborts the program whenever it gets to the default statement.
1111 * This is perhaps more difficult to read than necessary, and consequently
1113 * as the form above (with the minor difference that it also aborts the
1114 * program in release mode). It is written in all-caps because that makes
1115 * it stand out visually (and also because it is not actually a function,
1119 * template <
int dim>
1120 *
void ElasticProblem<dim>::output_results(const
unsigned int cycle) const
1125 * std::vector<std::string> solution_names;
1129 * solution_names.emplace_back(
"displacement");
1132 * solution_names.emplace_back(
"x_displacement");
1133 * solution_names.emplace_back(
"y_displacement");
1136 * solution_names.emplace_back(
"x_displacement");
1137 * solution_names.emplace_back(
"y_displacement");
1138 * solution_names.emplace_back(
"z_displacement");
1146 * After setting up the names
for the different components of the
1147 * solution vector, we can add the solution vector to the list of
1148 * data vectors scheduled
for output. Note that the following
1149 * function takes a vector of strings as
second argument, whereas
1150 * the one which we have used in all previous examples accepted a
1151 *
string there. (In fact, the function we had used before would
1152 * convert the single
string into a vector with only one element
1153 * and forwards that to the other function.)
1156 * data_out.add_data_vector(solution, solution_names);
1157 * data_out.build_patches();
1159 * std::ofstream output(
"solution-" + std::to_string(cycle) +
".vtk");
1160 * data_out.write_vtk(output);
1168 * <a name=
"step_8-ElasticProblemrun"></a>
1169 * <h4>ElasticProblem::run</h4>
1173 * The <code>
run</code> function does the same things as in @ref step_6
"step-6",
for
1174 * example. This time, we use the square [-1,1]^
d as domain, and we
refine
1175 * it globally four times before starting the
first iteration.
1179 * The reason
for refining is a bit accidental: we use the
QGauss
1180 * quadrature formula with two points in each direction
for integration of the
1181 * right hand side; that means that there are four quadrature points on each
1182 * cell (in 2d). If we only
refine the
initial grid once globally, then there
1183 * will be only four quadrature points in each direction on the
1184 * domain. However, the right hand side function was chosen to be rather
1185 * localized and in that
case, by pure chance, it happens that all quadrature
1186 * points lie at points where the right hand side function is zero (in
1187 * mathematical terms, the quadrature points happen to be at points outside
1188 * the <i>support</i> of the right hand side function). The right hand side
1189 * vector computed with quadrature will then contain only zeroes (even though
1190 * it would of course be nonzero
if we had computed the right hand side vector
1191 * exactly
using the integral) and the solution of the system of
1192 * equations is the zero vector, i.e., a finite element function that is zero
1193 * everywhere. In a sense, we
1194 * should not be surprised that
this is happening since we have chosen
1195 * an
initial grid that is totally unsuitable
for the problem at hand.
1199 * The unfortunate thing is that
if the discrete solution is
constant, then
1201 *
for each cell as well, and the call to
1202 * Triangulation::refine_and_coarsen_fixed_number() will not flag any cells
1203 * for refinement (why should it if the indicated error is zero for each
1204 * cell?). The grid in the next iteration will therefore consist of four
1205 * cells only as well, and the same problem occurs again.
1209 * The conclusion needs to be: while of course we will not choose the
1210 * initial grid to be well-suited for the accurate solution of the problem,
1211 * we must at least choose it such that it has the chance to capture the
1212 * important features of the solution. In this case, it needs to be able to
1213 * see the right hand side. Thus, we refine globally four times. (Any larger
1214 * number of global refinement steps would of course also work.)
1217 * template <
int dim>
1218 *
void ElasticProblem<dim>::run()
1220 *
for (
unsigned int cycle = 0; cycle < 8; ++cycle)
1222 * std::cout <<
"Cycle " << cycle <<
':' << std::endl;
1232 * std::cout <<
" Number of active cells: "
1237 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
1240 * assemble_system();
1242 * output_results(cycle);
1250 * <a name=
"step_8-Thecodemaincodefunction"></a>
1251 * <h3>The <code>main</code> function</h3>
1255 * After closing the <code>Step8</code>
namespace in the last line above, the
1256 * following is the main function of the program and is again exactly like in
1257 * @ref step_6
"step-6" (apart from the changed
class names, of course).
1264 * Step8::ElasticProblem<2> elastic_problem_2d;
1265 * elastic_problem_2d.run();
1267 *
catch (std::exception &exc)
1269 * std::cerr << std::endl
1271 * <<
"----------------------------------------------------"
1273 * std::cerr <<
"Exception on processing: " << std::endl
1274 * << exc.what() << std::endl
1275 * <<
"Aborting!" << std::endl
1276 * <<
"----------------------------------------------------"
1283 * std::cerr << std::endl
1285 * <<
"----------------------------------------------------"
1287 * std::cerr <<
"Unknown exception!" << std::endl
1288 * <<
"Aborting!" << std::endl
1289 * <<
"----------------------------------------------------"
1297<a name=
"step_8-Results"></a><h1>Results</h1>
1301There is not much to be said about the results of
this program, other than
1302that they look nice. All images were made
using VisIt from the
1303output files that the program wrote to disk. The
first two pictures show
1304the @f$x@f$- and @f$y@f$-displacements as
scalar components:
1309<img src=
"https://www.dealii.org/images/steps/developer/step-8.x.png" alt=
"">
1312<img src=
"https://www.dealii.org/images/steps/developer/step-8.y.png" alt=
"">
1318You can clearly see the sources of @f$x@f$-displacement around @f$x=0.5@f$ and
1319@f$x=-0.5@f$, and of @f$y@f$-displacement at the origin.
1321What one frequently would like to
do is to show the displacement as a vector
1322field, i.e., vectors that
for each point illustrate the direction and magnitude
1323of displacement. Unfortunately, that
's a bit more involved. To understand why
1324this is so, remember that we have just defined our finite element as a
1325collection of two components (in <code>dim=2</code> dimensions). Nowhere have
1326we said that this is not just a pressure and a concentration (two scalar
1327quantities) but that the two components actually are the parts of a
1328vector-valued quantity, namely the displacement. Absent this knowledge, the
1329DataOut class assumes that all individual variables we print are separate
1330scalars, and VisIt and Paraview then faithfully assume that this is indeed what it is. In
1331other words, once we have written the data as scalars, there is nothing in
1332these programs that allows us to paste these two scalar fields back together as a
1333vector field. Where we would have to attack this problem is at the root,
1334namely in <code>ElasticProblem::output_results()</code>. We won't
do so here but
1335instead refer the reader to the @ref step_22
"step-22" program where we show how to
do this
1336for a more
general situation. That said, we couldn
't help generating the data
1337anyway that would show how this would look if implemented as discussed in
1338@ref step_22 "step-22". The vector field then looks like this (VisIt and Paraview
1339randomly select a few
1340hundred vertices from which to draw the vectors; drawing them from each
1341individual vertex would make the picture unreadable):
1343<img src="https://www.dealii.org/images/steps/developer/step-8.vectors.png" alt="">
1346We note that one may have intuitively expected the
1347solution to be symmetric about the @f$x@f$- and @f$y@f$-axes since the @f$x@f$- and
1348@f$y@f$-forces are symmetric with respect to these axes. However, the force
1349considered as a vector is not symmetric and consequently neither is
1353<a name="step_8-PlainProg"></a>
1354<h1> The plain program</h1>
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
#define Assert(cond, exc)
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())
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
#define DEAL_II_NOT_IMPLEMENTED()
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
void refine(Triangulation< dim, spacedim > &tria, const Vector< Number > &criteria, const double threshold, const unsigned int max_to_mark=numbers::invalid_unsigned_int)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ general
No special properties.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
VectorType::value_type * begin(VectorType &V)
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)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation