404 *
double return_value = 0.0;
405 *
for (
unsigned int i = 0; i < dim; ++i)
406 * return_value += 4.0 *
std::pow(p(i), 4.0);
408 *
return return_value;
414 * As boundary
values, we choose @f$x^2+y^2@f$ in 2D, and @f$x^2+y^2+z^2@f$ in 3D. This
415 * happens to be
equal to the square of the vector from the origin to the
416 *
point at which we would like to evaluate the function, irrespective of the
417 * dimension. So that is what we
return:
421 *
double BoundaryValues<dim>::value(
const Point<dim> &p,
422 *
const unsigned int )
const
432 * <a name=
"ImplementationofthecodeStep4codeclass"></a>
433 * <h3>Implementation of the <code>Step4</code>
class</h3>
437 * Next
for the implementation of the
class template that makes use of the
438 *
functions above. As before, we will write everything as templates that have
439 * a formal parameter <code>dim</code> that we assume unknown at the time we
440 * define the
template functions. Only later, the compiler will find a
441 * declaration of <code>Step4@<2@></code> (in the <code>main</code> function,
442 * actually) and compile the entire
class with <code>dim</code> replaced by 2,
443 * a process referred to as `instantiation of a
template'. When doing so, it
444 * will also replace instances of <code>RightHandSide@<dim@></code> by
445 * <code>RightHandSide@<2@></code> and instantiate the latter class from the
450 * In fact, the compiler will also find a declaration <code>Step4@<3@></code>
451 * in <code>main()</code>. This will cause it to again go back to the general
452 * <code>Step4@<dim@></code> template, replace all occurrences of
453 * <code>dim</code>, this time by 3, and compile the class a second time. Note
454 * that the two instantiations <code>Step4@<2@></code> and
455 * <code>Step4@<3@></code> are completely independent classes; their only
456 * common feature is that they are both instantiated from the same general
457 * template, but they are not convertible into each other, for example, and
458 * share no code (both instantiations are compiled completely independently).
466 * <a name="Step4Step4"></a>
467 * <h4>Step4::Step4</h4>
471 * After this introduction, here is the constructor of the <code>Step4</code>
472 * class. It specifies the desired polynomial degree of the finite elements
473 * and associates the DoFHandler to the triangulation just as in the previous
474 * example program, @ref step_3 "step-3":
478 * Step4<dim>::Step4()
480 * , dof_handler(triangulation)
487 * <a name="Step4make_grid"></a>
488 * <h4>Step4::make_grid</h4>
492 * Grid creation is something inherently dimension dependent. However, as long
493 * as the domains are sufficiently similar in 2D or 3D, the library can
494 * abstract for you. In our case, we would like to again solve on the square
495 * @f$[-1,1]\times [-1,1]@f$ in 2D, or on the cube @f$[-1,1] \times [-1,1] \times
496 * [-1,1]@f$ in 3D; both can be termed GridGenerator::hyper_cube(), so we may
497 * use the same function in whatever dimension we are. Of course, the
498 * functions that create a hypercube in two and three dimensions are very much
499 * different, but that is something you need not care about. Let the library
500 * handle the difficult things.
504 * void Step4<dim>::make_grid()
506 * GridGenerator::hyper_cube(triangulation, -1, 1);
507 * triangulation.refine_global(4);
509 * std::cout << " Number of active cells: " << triangulation.n_active_cells()
511 * << " Total number of cells: " << triangulation.n_cells()
518 * <a name="Step4setup_system"></a>
519 * <h4>Step4::setup_system</h4>
523 * This function looks exactly like in the previous example, although it
524 * performs actions that in their details are quite different if
525 * <code>dim</code> happens to be 3. The only significant difference from a
526 * user's perspective is the number of cells resulting, which is much higher
527 * in three than in two space dimensions!
531 *
void Step4<dim>::setup_system()
533 * dof_handler.distribute_dofs(fe);
535 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
540 * sparsity_pattern.copy_from(dsp);
542 * system_matrix.reinit(sparsity_pattern);
544 * solution.reinit(dof_handler.n_dofs());
545 * system_rhs.reinit(dof_handler.n_dofs());
552 * <a name=
"Step4assemble_system"></a>
553 * <h4>Step4::assemble_system</h4>
557 * Unlike in the previous example, we would now like to use a non-
constant
558 * right hand side function and non-zero boundary
values. Both are tasks that
559 * are readily achieved with only a few
new lines of code in the assemblage of
560 * the
matrix and right hand side.
564 * More interesting, though, is the way we
assemble matrix and right hand side
565 * vector dimension independently: there is simply no difference to the
566 * two-dimensional
case. Since the important objects used in
this function
567 * (quadrature formula,
FEValues) depend on the dimension by way of a
template
568 * parameter as well, they can take care of setting up properly everything
for
569 * the dimension
for which
this function is compiled. By declaring all classes
570 * which might depend on the dimension
using a
template parameter, the library
571 * can make nearly all work
for you and you don
't have to care about most
576 * void Step4<dim>::assemble_system()
578 * QGauss<dim> quadrature_formula(fe.degree + 1);
582 * We wanted to have a non-constant right hand side, so we use an object of
583 * the class declared above to generate the necessary data. Since this right
584 * hand side object is only used locally in the present function, we declare
585 * it here as a local variable:
588 * RightHandSide<dim> right_hand_side;
592 * Compared to the previous example, in order to evaluate the non-constant
593 * right hand side function we now also need the quadrature points on the
594 * cell we are presently on (previously, we only required values and
595 * gradients of the shape function from the FEValues object, as well as the
596 * quadrature weights, FEValues::JxW() ). We can tell the FEValues object to
597 * do for us by also giving it the #update_quadrature_points flag:
600 * FEValues<dim> fe_values(fe,
601 * quadrature_formula,
602 * update_values | update_gradients |
603 * update_quadrature_points | update_JxW_values);
607 * We then again define the same abbreviation as in the previous program.
608 * The value of this variable of course depends on the dimension which we
609 * are presently using, but the FiniteElement class does all the necessary
610 * work for you and you don't have to care about the dimension dependent
614 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
619 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
623 * Next, we again have to
loop over all cells and
assemble local
624 * contributions. Note, that a cell is a quadrilateral in two space
625 * dimensions, but a hexahedron in 3D. In fact, the
626 * <code>active_cell_iterator</code> data type is something different,
627 * depending on the dimension we are in, but to the
outside world they look
628 * alike and you will probably never see a difference. In any
case, the real
629 * type is hidden by
using `
auto`:
632 *
for (
const auto &cell : dof_handler.active_cell_iterators())
634 * fe_values.reinit(cell);
640 * Now we have to
assemble the local
matrix and right hand side. This is
641 * done exactly like in the previous example, but now we revert the
642 * order of the loops (which we can safely
do since they are independent
643 * of each other) and
merge the loops
for the local
matrix and the local
644 * vector as far as possible to make things a bit faster.
648 * Assembling the right hand side presents the only significant
649 * difference to how we did things in @ref step_3
"step-3": Instead of
using a
650 *
constant right hand side with
value 1, we use the
object representing
651 * the right hand side and evaluate it at the quadrature points:
654 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
655 *
for (
const unsigned int i : fe_values.dof_indices())
657 *
for (
const unsigned int j : fe_values.dof_indices())
659 * (fe_values.shape_grad(i, q_index) *
660 * fe_values.shape_grad(j, q_index) *
661 * fe_values.JxW(q_index));
663 *
const auto &x_q = fe_values.quadrature_point(q_index);
664 * cell_rhs(i) += (fe_values.shape_value(i, q_index) *
665 * right_hand_side.value(x_q) *
666 * fe_values.JxW(q_index));
670 * As a
final remark to these loops: when we
assemble the local
671 * contributions into <code>
cell_matrix(i,j)</code>, we have to multiply
674 * multiply it with the
scalar weights JxW. This is what actually
675 * happens: <code>fe_values.shape_grad(i,q_index)</code> returns a
676 * <code>dim</code> dimensional vector, represented by a
677 * <code>
Tensor@<1,dim@></code> object, and the
operator* that
678 * multiplies it with the result of
679 * <code>fe_values.shape_grad(j,q_index)</code> makes sure that the
680 * <code>dim</code> components of the two vectors are properly
681 * contracted, and the result is a
scalar floating
point number that
682 * then is multiplied with the weights. Internally,
this operator* makes
683 * sure that
this happens correctly
for all <code>dim</code> components
684 * of the vectors, whether <code>dim</code> be 2, 3, or any other space
685 * dimension; from a user
's perspective, this is not something worth
686 * bothering with, however, making things a lot simpler if one wants to
687 * write code dimension independently.
691 * With the local systems assembled, the transfer into the global matrix
692 * and right hand side is done exactly as before, but here we have again
693 * merged some loops for efficiency:
696 * cell->get_dof_indices(local_dof_indices);
697 * for (const unsigned int i : fe_values.dof_indices())
699 * for (const unsigned int j : fe_values.dof_indices())
700 * system_matrix.add(local_dof_indices[i],
701 * local_dof_indices[j],
702 * cell_matrix(i, j));
704 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
710 * As the final step in this function, we wanted to have non-homogeneous
711 * boundary values in this example, unlike the one before. This is a simple
712 * task, we only have to replace the Functions::ZeroFunction used there by an
713 * object of the class which describes the boundary values we would like to
714 * use (i.e. the <code>BoundaryValues</code> class declared above):
718 * The function VectorTools::interpolate_boundary_values() will only work
719 * on faces that have been marked with boundary indicator 0 (because that's
720 * what we say the function should work on with the
second argument below).
721 * If there are faces with boundary
id other than 0, then the function
722 * interpolate_boundary_values will
do nothing on these faces. For
723 * the Laplace equation doing nothing is equivalent to assuming that
724 * on those parts of the boundary a zero Neumann boundary condition holds.
727 * std::map<types::global_dof_index, double> boundary_values;
730 * BoundaryValues<dim>(),
742 * <a name=
"Step4solve"></a>
743 * <h4>Step4::solve</h4>
747 * Solving the linear system of equations is something that looks almost
748 * identical in most programs. In particular, it is dimension independent, so
749 *
this function is copied verbatim from the previous example.
753 *
void Step4<dim>::solve()
761 * We have made one addition, though: since we suppress output from the
762 * linear solvers, we have to print the number of iterations by hand.
765 * std::cout <<
" " << solver_control.last_step()
766 * <<
" CG iterations needed to obtain convergence." << std::endl;
773 * <a name=
"Step4output_results"></a>
774 * <h4>Step4::output_results</h4>
778 * This function also does what the respective one did in @ref step_3
"step-3". No changes
779 * here
for dimension independence either.
783 * Since the program will
run both 2
d and 3
d versions of the Laplace solver,
784 * we use the dimension in the filename to generate distinct filenames
for
785 * each
run (in a better program, one would check whether <code>dim</code> can
786 * have other values than 2 or 3, but we neglect
this here
for the sake of
791 *
void Step4<dim>::output_results() const
800 * std::ofstream output(dim == 2 ?
"solution-2d.vtk" :
"solution-3d.vtk");
809 * <a name=
"Step4run"></a>
810 * <h4>Step4::run</h4>
814 * This is the function which has the top-
level control over everything. Apart
815 * from one line of additional output, it is the same as
for the previous
820 *
void Step4<dim>::run()
822 * std::cout <<
"Solving problem in " << dim <<
" space dimensions."
836 * <a name=
"Thecodemaincodefunction"></a>
837 * <h3>The <code>main</code> function</h3>
841 * And
this is the main function. It also looks mostly like in @ref step_3
"step-3", but
if
842 * you look at the code below, note how we
first create a variable of type
843 * <code>Step4@<2@></code> (forcing the compiler to compile the
class template
844 * with <code>dim</code> replaced by <code>2</code>) and run a 2d simulation,
845 * and then we
do the whole thing over in 3
d.
849 * In practice,
this is probably not what you would
do very frequently (you
850 * probably either want to solve a 2d problem, or one in 3d, but not both at
851 * the same time). However, it demonstrates the mechanism by which we can
852 * simply change which dimension we want in a single place, and thereby force
853 * the compiler to recompile the dimension independent
class templates for the
854 * dimension we request. The emphasis here lies on the fact that we only need
855 * to change a single place. This makes it rather trivial to debug the program
856 * in 2
d where computations are fast, and then
switch a single place to a 3 to
857 *
run the much more computing intensive program in 3
d for `real
'
862 * Each of the two blocks is enclosed in braces to make sure that the
863 * <code>laplace_problem_2d</code> variable goes out of scope (and releases
864 * the memory it holds) before we move on to allocate memory for the 3d
865 * case. Without the additional braces, the <code>laplace_problem_2d</code>
866 * variable would only be destroyed at the end of the function, i.e. after
867 * running the 3d problem, and would needlessly hog memory while the 3d run
868 * could actually use it.
874 * Step4<2> laplace_problem_2d;
875 * laplace_problem_2d.run();
879 * Step4<3> laplace_problem_3d;
880 * laplace_problem_3d.run();
886<a name="Results"></a><h1>Results</h1>
890The output of the program looks as follows (the number of iterations
891may vary by one or two, depending on your computer, since this is
892often dependent on the round-off accuracy of floating point
893operations, which differs between processors):
895Solving problem in 2 space dimensions.
896 Number of active cells: 256
897 Total number of cells: 341
898 Number of degrees of freedom: 289
899 26 CG iterations needed to obtain convergence.
900Solving problem in 3 space dimensions.
901 Number of active cells: 4096
902 Total number of cells: 4681
903 Number of degrees of freedom: 4913
904 30 CG iterations needed to obtain convergence.
906It is obvious that in three spatial dimensions the number of cells and
907therefore also the number of degrees of freedom is
908much higher. What cannot be seen here, is that besides this higher
909number of rows and columns in the matrix, there are also significantly
910more entries per row of the matrix in three space
911dimensions. Together, this leads to a much higher numerical effort for
912solving the system of equation, which you can feel in the run time of the two
913solution steps when you actually run the program.
917The program produces two files: <code>solution-2d.vtk</code> and
918<code>solution-3d.vtk</code>, which can be viewed using the programs
919VisIt or Paraview (in case you do not have these programs, you can easily
921output format in the program to something which you can view more
922easily). Visualizing solutions is a bit of an art, but it can also be fun, so
923you should play around with your favorite visualization tool to get familiar
924with its functionality. Here's what I have come up with
for the 2
d solution:
927 <img src=
"https://www.dealii.org/images/steps/developer/step-4.solution-2d.png" alt=
"">
930(See also <a href=
"http://www.math.colostate.edu/~bangerth/videos.676.11.html">video lecture 11</a>, <a href=
"http://www.math.colostate.edu/~bangerth/videos.676.32.html">video lecture 32</a>.)
931The picture shows the solution of the problem under consideration as
932a 3D plot. As can be seen, the solution is almost flat in the interior
933of the domain and has a higher curvature near the boundary. This, of
934course, is due to the fact that for Laplace's equation the curvature
935of the solution is
equal to the right hand side and that was chosen as
936a quartic polynomial which is nearly zero in the interior and is only
937rising sharply when approaching the boundaries of the domain; the
938maximal values of the right hand side function are at the corners of
939the domain, where also the solution is moving most rapidly.
940It is also nice to see that the solution follows the desired quadratic
941boundary values along the boundaries of the domain.
942It can also be useful to verify a computed solution against an analytical
943solution. For an explanation of this technique, see @ref step_7 "step-7".
945On the other hand, even though the picture does not show the mesh lines
946explicitly, you can see them as little kinks in the solution. This clearly
947indicates that the solution hasn't been computed to very high accuracy and
948that to get a better solution, we may have to compute on a finer mesh.
950In three spatial dimensions, visualization is a bit more difficult. The left
951picture shows the solution and the mesh it was computed on on the surface of
952the domain. This is nice, but it has the drawback that it completely hides
953what is happening on the inside. The picture on the right is an attempt at
954visualizing the interior as well, by showing surfaces where the solution has
955constant values (as indicated by the legend at the top left). Isosurface
956pictures look best if one makes the individual surfaces slightly transparent
957so that it is possible to see through them and see what's behind.
959<table width="60%" align="
center">
971A final remark on visualization: the idea of visualization is to give insight,
972which is not the same as displaying information. In particular, it is easy to
973overload a picture with information, but while it shows more information it
974makes it also more difficult to glean insight. As an example, the program I
975used to generate these pictures, VisIt, by default puts tick marks on every
976axis, puts a big fat label "X Axis" on the @f$x@f$ axis and similar for the other
977axes, shows the file name from which the data was taken in the top left and
978the name of the user doing so and the time and date on the bottom right. None
980here: the axes are equally easy to make out because the tripod at the bottom
981left is still visible, and we know from the program that the domain is
982@f$[-1,1]^3@f$, so there is no need for tick marks. As a consequence, I have
983switched off all the extraneous stuff in the picture: the art of visualization
984is to reduce the picture to those parts that are important to see what one
985wants to see, but no more.
989<a name="extensions"></a>
990<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
994Essentially the possibilities for playing around with the program are the same
995as for the previous one, except that they will now also apply to the 3d
996case. For inspiration read up on <a href="step_3.html
#extensions"
997target="body">possible extensions in the documentation of step 3</a>.
1000<a name="PlainProg"></a>
1001<h1> The plain program</h1>
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void add_data_vector(const VectorType &data, const std::vector< std::string > &names, const DataVectorType type=type_automatic, const std::vector< DataComponentInterpretation::DataComponentInterpretation > &data_component_interpretation={})
virtual void build_patches(const unsigned int n_subdivisions=0)
numbers::NumberTraits< Number >::real_type square() const
void write_vtk(std::ostream &out) const
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternType &sparsity_pattern, const AffineConstraints< number > &constraints=AffineConstraints< number >(), const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ matrix
Contents is actually a matrix.
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)
void run(const Iterator &begin, const typename identity< Iterator >::type &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)
constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)