409 * <a name=
"step_5-Step5setup_system"></a>
410 * <h4>Step5::setup_system</h4>
414 * This is the function <code>make_grid</code> from the previous
415 * example, minus the generation of the grid. Everything
else is unchanged:
419 *
void Step5<dim>::setup_system()
421 * dof_handler.distribute_dofs(fe);
423 * std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
428 * sparsity_pattern.copy_from(dsp);
430 * system_matrix.reinit(sparsity_pattern);
432 * solution.reinit(dof_handler.n_dofs());
433 * system_rhs.reinit(dof_handler.n_dofs());
441 * <a name=
"step_5-Step5assemble_system"></a>
442 * <h4>Step5::assemble_system</h4>
446 * As in the previous examples,
this function is not changed much with regard
447 * to its functionality, but there are still some optimizations which we will
448 * show. For
this, it is important to note that
if efficient solvers are used
449 * (such as the preconditioned CG method), assembling the matrix and right hand
450 * side can take a comparable time, and you should think about
using one or
451 * two optimizations at some places.
455 * The
first parts of the function are completely unchanged from before:
459 * void Step5<dim>::assemble_system()
461 * const
QGauss<dim> quadrature_formula(fe.degree + 1);
464 * quadrature_formula,
468 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
473 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
477 * Next is the typical
loop over all cells to compute local contributions
478 * and then to transfer them into the global
matrix and vector. The only
479 * change in
this part, compared to @ref step_4
"step-4", is that we will use the
480 * <code>coefficient()</code> function defined above to compute the
484 *
for (
const auto &cell : dof_handler.active_cell_iterators())
491 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
493 * const double current_coefficient =
494 * coefficient(fe_values.quadrature_point(q_index));
495 *
for (
const unsigned int i : fe_values.dof_indices())
497 * for (const unsigned
int j : fe_values.dof_indices())
499 * (current_coefficient *
500 * fe_values.shape_grad(i, q_index) *
501 * fe_values.shape_grad(j, q_index) *
502 * fe_values.JxW(q_index));
504 * cell_rhs(i) += (fe_values.shape_value(i, q_index) *
506 * fe_values.JxW(q_index));
511 * cell->get_dof_indices(local_dof_indices);
512 *
for (
const unsigned int i : fe_values.dof_indices())
514 * for (const unsigned
int j : fe_values.dof_indices())
515 * system_matrix.add(local_dof_indices[i],
516 * local_dof_indices[j],
519 * system_rhs(local_dof_indices[i]) += cell_rhs(i);
525 * With the
matrix so built, we use zero boundary
values again:
528 * std::map<types::global_dof_index, double> boundary_values;
543 * <a name=
"step_5-Step5solve"></a>
544 * <h4>Step5::solve</h4>
548 * The solution process again looks mostly like in the previous
549 * examples. However, we will now use a preconditioned conjugate
gradient
550 * algorithm. It is not very difficult to make
this change. In fact, the only
551 * thing we have to alter is that we need an
object which will act as a
552 * preconditioner. We will use SSOR (symmetric successive overrelaxation),
553 * with a relaxation factor of 1.2. For
this purpose, the
554 * <code>
SparseMatrix</code>
class has a function which does one SSOR step,
555 * and we need to
package the address of this function together with the
556 * matrix on which it should act (which is the matrix to be inverted) and the
557 * relaxation factor into one object. The <code>PreconditionSSOR</code> class
558 * does this for us. (<code>PreconditionSSOR</code> class takes a template
559 * argument denoting the matrix type it is supposed to work on. The default
560 * value is <code>SparseMatrix@<double@></code>, which is exactly what we need
561 * here, so we simply stick with the default and do not specify anything in
562 * the angle brackets.)
566 * Note that for the present case, SSOR doesn't really perform much better
567 * than most other preconditioners (though better than no preconditioning at
568 * all). A brief comparison of different preconditioners is presented in the
569 * Results section of the next tutorial program, @ref step_6 "step-6".
573 * With this, the rest of the function is trivial: instead of the
574 * <code>PreconditionIdentity</code> object we have created before, we now use
575 * the preconditioner we have declared, and the CG solver will do the rest for
580 * void Step5<dim>::solve()
582 * SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
586 * preconditioner.
initialize(system_matrix, 1.2);
588 * solver.solve(system_matrix, solution, system_rhs, preconditioner);
590 * std::cout <<
" " << solver_control.last_step()
591 * <<
" CG iterations needed to obtain convergence." << std::endl;
598 * <a name=
"step_5-Step5output_resultsandsettingoutputflags"></a>
599 * <h4>Step5::output_results and setting output flags</h4>
603 * Writing output to a file is mostly the same as
for the previous tutorial.
604 * The only difference is that we now need to construct a different filename
605 *
for each refinement cycle.
609 * The function writes the output in VTU format, a variation of the VTK format
610 * that
requires less disk space because it compresses the data. Of course,
611 * there are many other formats supported by the
DataOut class if you
612 * desire to use a program
for visualization that doesn
't understand
617 * void Step5<dim>::output_results(const unsigned int cycle) const
619 * DataOut<dim> data_out;
621 * data_out.attach_dof_handler(dof_handler);
622 * data_out.add_data_vector(solution, "solution");
624 * data_out.build_patches();
626 * std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
627 * data_out.write_vtu(output);
635 * <a name="step_5-Step5run"></a>
636 * <h4>Step5::run</h4>
640 * The second to last thing in this program is the definition of the
641 * <code>run()</code> function. In contrast to the previous programs, we will
642 * compute on a sequence of meshes that after each iteration is globally
643 * refined. The function therefore consists of a loop over 6 cycles. In each
644 * cycle, we first print the cycle number, and then have to decide what to do
645 * with the mesh. If this is not the first cycle, we simply refine the
646 * existing mesh once globally. Before running through these cycles, however,
647 * we have to generate a mesh:
651 * In previous examples, we have already used some of the functions from the
652 * <code>GridGenerator</code> class. Here we would like to read a grid from a
653 * file where the cells are stored and which may originate from someone else,
654 * or may be the product of a mesh generator tool.
658 * In order to read a grid from a file, we generate an object of data type
659 * GridIn and associate the triangulation to it (i.e. we tell it to fill our
660 * triangulation object when we ask it to read the file). Then we open the
661 * respective file and initialize the triangulation with the data in the file :
665 * void Step5<dim>::run()
667 * GridIn<dim> grid_in;
668 * grid_in.attach_triangulation(triangulation);
669 * std::ifstream input_file("circle-grid.inp");
672 * We would now like to read the file. However, the input file is only for a
673 * two-dimensional triangulation, while this function is a template for
674 * arbitrary dimension. Since this is only a demonstration program, we will
675 * not use different input files for the different dimensions, but rather
676 * quickly kill the whole program if we are not in 2d. Of course, since the
677 * main function below assumes that we are working in two dimensions we
678 * could skip this check, in this version of the program, without any ill
683 * It turns out that perhaps 90 per cent of programming errors are invalid
684 * function parameters such as invalid array sizes, etc., so we use
685 * assertions heavily throughout deal.II to catch such mistakes. For this,
686 * the <code>Assert</code> macro is a good choice, since it makes sure that
687 * the condition which is given as first argument is valid, and if not
688 * throws an exception (its second argument) which will usually terminate
689 * the program giving information where the error occurred and what the
690 * reason was. (A longer discussion of what exactly the @p Assert macro
691 * does can be found in the @ref Exceptions "exception documentation topic".)
692 * This generally reduces the time to find programming errors
693 * dramatically and we have found assertions an invaluable means to program
698 * On the other hand, all these checks (there are over 10,000 of them in the
699 * library at present) should not slow down the program too much if you want
700 * to do large computations. To this end, the <code>Assert</code> macro is
701 * only used in debug mode and expands to nothing if in optimized
702 * mode. Therefore, while you test your program on small problems and debug
703 * it, the assertions will tell you where the problems are. Once your
704 * program is stable, you can switch off debugging and the program will run
705 * your real computations without the assertions and at maximum speed. More
706 * precisely: turning off all the checks in the library (which prevent you
707 * from calling functions with wrong arguments, walking off of arrays, etc.)
708 * by compiling your program in optimized mode usually makes things run
709 * about four times faster. Even though optimized programs are more
710 * performant, you should always develop in debug mode since it allows
711 * the library to find lots of common programming errors automatically. For
712 * those who want to try: The way to switch from debug mode to optimized
713 * mode is to recompile your program with the command <code>make
714 * release</code>. The output of the <code>make</code> program should now
715 * indicate to you that the program is now compiled in optimized mode, and
716 * it will later also be linked to libraries that have been compiled for
717 * optimized mode. In order to switch back to debug mode, simply recompile
718 * with the command <code>make debug</code>.
721 * Assert(dim == 2, ExcNotImplemented());
724 * ExcNotImplemented is a globally defined exception, which may be thrown
725 * whenever a piece of code has simply not been implemented for a case
726 * other than the condition checked in the assertion. Here, it would not
727 * be difficult to simply implement reading a *different* mesh file that
728 * contains a description of a 1d or 3d geometry, but this has not (yet)
729 * been implemented and so the exception is appropriate.
733 * Usually, one would like to use more specific
734 * exception classes, and particular in this case one would of course try
735 * to do something else if <code>dim</code> is not equal to two, e.g. create
736 * a grid using library functions. Aborting a program is usually not a good
737 * idea and assertions should really only be used for exceptional cases
738 * which should not occur, but might due to stupidity of the programmer,
739 * user, or someone else.
743 * So if we got past the assertion, we know that dim==2, and we can now
744 * actually read the grid. It is in UCD (unstructured cell data) format
745 * (though the convention is to use the suffix <code>inp</code> for UCD
749 * grid_in.read_ucd(input_file);
752 * If you like to use another input format, you have to use one of the other
753 * <code>grid_in.read_xxx</code> function. (See the documentation of the
754 * <code>GridIn</code> class to find out what input formats are presently
759 * The grid in the file describes a circle. Therefore we have to use a
760 * manifold object which tells the triangulation where to put new points on
761 * the boundary when the grid is refined. Unlike @ref step_1 "step-1", since GridIn does
762 * not know that the domain has a circular boundary (unlike
763 * GridGenerator::hyper_shell) we have to explicitly attach a manifold to
764 * the boundary after creating the triangulation to get the correct result
765 * when we refine the mesh.
768 * const SphericalManifold<dim> boundary;
769 * triangulation.set_all_manifold_ids_on_boundary(0);
770 * triangulation.set_manifold(0, boundary);
772 * for (unsigned int cycle = 0; cycle < 6; ++cycle)
774 * std::cout << "Cycle " << cycle << ':
' << std::endl;
777 * triangulation.refine_global(1);
781 * Now that we have a mesh for sure, we write some output and do all the
782 * things that we have already seen in the previous examples.
785 * std::cout << " Number of active cells: "
786 * << triangulation.n_active_cells()
788 * << " Total number of cells: "
789 * << triangulation.n_cells()
795 * output_results(cycle);
803 * <a name="step_5-Thecodemaincodefunction"></a>
804 * <h3>The <code>main</code> function</h3>
808 * The main function looks mostly like the one in the previous example, so we
809 * won't comment on it further:
814 * Step5<2> laplace_problem_2d;
815 * laplace_problem_2d.run();
819<a name=
"step_5-Results"></a><h1>Results</h1>
823Here is the console output:
826 Number of active cells: 20
827 Total number of cells: 20
828 Number of degrees of freedom: 25
829 8 CG iterations needed to obtain convergence.
831 Number of active cells: 80
832 Total number of cells: 100
833 Number of degrees of freedom: 89
834 12 CG iterations needed to obtain convergence.
836 Number of active cells: 320
837 Total number of cells: 420
838 Number of degrees of freedom: 337
839 21 CG iterations needed to obtain convergence.
841 Number of active cells: 1280
842 Total number of cells: 1700
843 Number of degrees of freedom: 1313
844 38 CG iterations needed to obtain convergence.
846 Number of active cells: 5120
847 Total number of cells: 6820
848 Number of degrees of freedom: 5185
849 70 CG iterations needed to obtain convergence.
851 Number of active cells: 20480
852 Total number of cells: 27300
853 Number of degrees of freedom: 20609
854 136 CG iterations needed to obtain convergence.
859In each cycle, the number of cells quadruples and the number of CG
860iterations roughly doubles.
861Also, in each cycle, the program writes one output graphic file in VTU
862format. They are depicted in the following:
867 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-0-r9.2.png" alt=
"">
870 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-1-r9.2.png" alt=
"">
875 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-2-r9.2.png" alt=
"">
878 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-3-r9.2.png" alt=
"">
883 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-4-r9.2.png" alt=
"">
886 <img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-5-r9.2.png" alt=
"">
893Due to the variable coefficient (the curvature there is reduced by the
894same factor by which the coefficient is increased), the top region of
895the solution is flattened. The
gradient of the solution is
896discontinuous along the interface, although
this is not very clearly
897visible in the pictures above. We will look at
this in more detail in
900The pictures also show that the solution computed by
this program is
901actually pretty wrong on a very coarse mesh (its magnitude is
902wrong). That
's because no numerical method guarantees that the solution
903on a coarse mesh is particularly accurate -- but we know that the
904solution <i>converges</i> to the exact solution, and indeed you can
905see how the solutions from one mesh to the next seem to not change
906very much any more at the end.
909<a name="step_5-PlainProg"></a>
910<h1> The plain program</h1>
void initialize(const MatrixType &A, const AdditionalData ¶meters=AdditionalData())
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)
@ update_values
Shape function values.
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
@ 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)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation