410 * <a name=
"step_5-Step5setup_system"></a>
422 *
dof_handler.distribute_dofs(fe);
424 *
std::cout <<
" Number of degrees of freedom: " << dof_handler.n_dofs()
429 *
sparsity_pattern.copy_from(
dsp);
431 *
system_matrix.reinit(sparsity_pattern);
433 *
solution.reinit(dof_handler.n_dofs());
442 * <a name=
"step_5-Step5assemble_system"></a>
469 *
const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
474 *
std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
485 *
for (
const auto &cell : dof_handler.active_cell_iterators())
492 *
for (
const unsigned int q_index : fe_values.quadrature_point_indices())
495 *
coefficient(fe_values.quadrature_point(q_index));
496 *
for (
const unsigned int i : fe_values.dof_indices())
501 *
fe_values.shape_grad(i, q_index) *
502 *
fe_values.shape_grad(
j, q_index) *
503 *
fe_values.JxW(q_index));
505 *
cell_rhs(i) += (fe_values.shape_value(i, q_index) *
507 *
fe_values.JxW(q_index));
512 *
cell->get_dof_indices(local_dof_indices);
513 *
for (
const unsigned int i : fe_values.dof_indices())
516 *
system_matrix.add(local_dof_indices[i],
517 *
local_dof_indices[
j],
544 * <a name=
"step_5-Step5solve"></a>
545 * <
h4>Step5::solve</
h4>
556 *
and we need to package the address of this function together with the
557 * matrix on which it should act (which is the matrix to be inverted) and the
558 * relaxation factor into one object. The <code>PreconditionSSOR</code> class
559 * does this for us. (<code>PreconditionSSOR</code> class takes a template
560 * argument denoting the matrix type it is supposed to work on. The default
561 * value is <code>SparseMatrix@<double@></code>, which is exactly what we need
562 * here, so we simply stick with the default and do not specify anything in
563 * the angle brackets.)
567 * Note that for the present case, SSOR doesn't really perform much better
568 * than most other preconditioners (though better than no preconditioning at
569 * all). A brief comparison of different preconditioners is presented in the
570 * Results section of the next tutorial program, @ref step_6 "step-6".
574 * With this, the rest of the function is trivial: instead of the
575 * <code>PreconditionIdentity</code> object we have created before, we now use
576 * the preconditioner we have declared, and the CG solver will do the rest for
581 * void Step5<dim>::solve()
583 * SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
587 *
preconditioner.initialize(system_matrix, 1.2);
589 *
solver.solve(system_matrix, solution,
system_rhs, preconditioner);
591 *
std::cout <<
" " << solver_control.last_step()
592 *
<<
" CG iterations needed to obtain convergence." << std::endl;
599 * <a name=
"step_5-Step5output_resultsandsettingoutputflags"></a>
618 * void Step5<dim>::output_results(const unsigned int cycle) const
620 * DataOut<dim> data_out;
622 * data_out.attach_dof_handler(dof_handler);
623 * data_out.add_data_vector(solution, "solution");
625 * data_out.build_patches();
627 * std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
628 * data_out.write_vtu(output);
636 * <a name="step_5-Step5run"></a>
637 * <h4>Step5::run</h4>
641 * The second to last thing in this program is the definition of the
642 * <code>run()</code> function. In contrast to the previous programs, we will
643 * compute on a sequence of meshes that after each iteration is globally
644 * refined. The function therefore consists of a loop over 6 cycles. In each
645 * cycle, we first print the cycle number, and then have to decide what to do
646 * with the mesh. If this is not the first cycle, we simply refine the
647 * existing mesh once globally. Before running through these cycles, however,
648 * we have to generate a mesh:
652 * In previous examples, we have already used some of the functions from the
653 * <code>GridGenerator</code> class. Here we would like to read a grid from a
654 * file where the cells are stored and which may originate from someone else,
655 * or may be the product of a mesh generator tool.
659 * In order to read a grid from a file, we generate an object of data type
660 * GridIn and associate the triangulation to it (i.e. we tell it to fill our
661 * triangulation object when we ask it to read the file). Then we open the
662 * respective file and initialize the triangulation with the data in the file :
666 * void Step5<dim>::run()
668 * GridIn<dim> grid_in;
669 * grid_in.attach_triangulation(triangulation);
670 * std::ifstream input_file("circle-grid.inp");
673 * We would now like to read the file. However, the input file is only for a
674 * two-dimensional triangulation, while this function is a template for
675 * arbitrary dimension. Since this is only a demonstration program, we will
676 * not use different input files for the different dimensions, but rather
677 * quickly kill the whole program if we are not in 2d. Of course, since the
678 * main function below assumes that we are working in two dimensions we
679 * could skip this check, in this version of the program, without any ill
684 * It turns out that perhaps 90 per cent of programming errors are invalid
685 * function parameters such as invalid array sizes, etc., so we use
686 * assertions heavily throughout deal.II to catch such mistakes. For this,
687 * the <code>Assert</code> macro is a good choice, since it makes sure that
688 * the condition which is given as first argument is valid, and if not
689 * throws an exception (its second argument) which will usually terminate
690 * the program giving information where the error occurred and what the
691 * reason was. (A longer discussion of what exactly the @p Assert macro
692 * does can be found in the @ref Exceptions "exception documentation topic".)
693 * This generally reduces the time to find programming errors
694 * dramatically and we have found assertions an invaluable means to program
699 * On the other hand, all these checks (there are over 10,000 of them in the
700 * library at present) should not slow down the program too much if you want
701 * to do large computations. To this end, the <code>Assert</code> macro is
702 * only used in debug mode and expands to nothing if in optimized
703 * mode. Therefore, while you test your program on small problems and debug
704 * it, the assertions will tell you where the problems are. Once your
705 * program is stable, you can switch off debugging and the program will run
706 * your real computations without the assertions and at maximum speed. More
707 * precisely: turning off all the checks in the library (which prevent you
708 * from calling functions with wrong arguments, walking off of arrays, etc.)
709 * by compiling your program in optimized mode usually makes things run
710 * about four times faster. Even though optimized programs are more
711 * performant, you should always develop in debug mode since it allows
712 * the library to find lots of common programming errors automatically. For
713 * those who want to try: The way to switch from debug mode to optimized
714 * mode is to recompile your program with the command <code>make
715 * release</code>. The output of the <code>make</code> program should now
716 * indicate to you that the program is now compiled in optimized mode, and
717 * it will later also be linked to libraries that have been compiled for
718 * optimized mode. In order to switch back to debug mode, simply recompile
719 * with the command <code>make debug</code>.
722 * Assert(dim == 2, ExcNotImplemented());
725 * ExcNotImplemented is a globally defined exception, which may be thrown
726 * whenever a piece of code has simply not been implemented for a case
727 * other than the condition checked in the assertion. Here, it would not
728 * be difficult to simply implement reading a *different* mesh file that
729 * contains a description of a 1d or 3d geometry, but this has not (yet)
730 * been implemented and so the exception is appropriate.
734 * Usually, one would like to use more specific
735 * exception classes, and particular in this case one would of course try
736 * to do something else if <code>dim</code> is not equal to two, e.g. create
737 * a grid using library functions. Aborting a program is usually not a good
738 * idea and assertions should really only be used for exceptional cases
739 * which should not occur, but might due to stupidity of the programmer,
740 * user, or someone else.
744 * So if we got past the assertion, we know that dim==2, and we can now
745 * actually read the grid. It is in UCD (unstructured cell data) format
746 * (though the convention is to use the suffix <code>inp</code> for UCD
750 * grid_in.read_ucd(input_file);
753 * If you like to use another input format, you have to use one of the other
754 * <code>grid_in.read_xxx</code> function. (See the documentation of the
755 * <code>GridIn</code> class to find out what input formats are presently
760 * The grid in the file describes a circle. Therefore we have to use a
761 * manifold object which tells the triangulation where to put new points on
762 * the boundary when the grid is refined. Unlike @ref step_1 "step-1", since GridIn does
763 * not know that the domain has a circular boundary (unlike
764 * GridGenerator::hyper_shell) we have to explicitly attach a manifold to
765 * the boundary after creating the triangulation to get the correct result
766 * when we refine the mesh.
769 * const SphericalManifold<dim> boundary;
770 * triangulation.set_all_manifold_ids_on_boundary(0);
771 * triangulation.set_manifold(0, boundary);
773 * for (unsigned int cycle = 0; cycle < 6; ++cycle)
775 * std::cout << "Cycle " << cycle << ':
' << std::endl;
778 * triangulation.refine_global(1);
782 * Now that we have a mesh for sure, we write some output and do all the
783 * things that we have already seen in the previous examples.
786 * std::cout << " Number of active cells: "
787 * << triangulation.n_active_cells()
789 * << " Total number of cells: "
790 * << triangulation.n_cells()
796 * output_results(cycle);
804 * <a name="step_5-Thecodemaincodefunction"></a>
805 * <h3>The <code>main</code> function</h3>
809 * The main function looks mostly like the one in the previous example, so we
868 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-0-r9.2.png" alt=
"">
871 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-1-r9.2.png" alt=
"">
876 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-2-r9.2.png" alt=
"">
879 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-3-r9.2.png" alt=
"">
884 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-4-r9.2.png" alt=
"">
887 <
img src=
"https://www.dealii.org/images/steps/developer/step-5.solution-5-r9.2.png" alt=
"">
903wrong).
That's because no numerical method guarantees that the solution
904on a coarse mesh is particularly accurate -- but we know that the
905solution <i>converges</i> to the exact solution, and indeed you can
906see how the solutions from one mesh to the next seem to not change
907very much any more at the end.
910<a name="step_5-PlainProg"></a>
911<h1> The plain program</h1>
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.
std::vector< index_type > data
@ matrix
Contents is actually a matrix.
constexpr types::blas_int zero
constexpr types::blas_int one
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