Reference documentation for deal.II version 9.6.0
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
step-5.h
Go to the documentation of this file.
1
399 1)
400 *   , dof_handler(triangulation)
401 *   {}
402 *  
403 *  
404 *  
405 * @endcode
406 *
407 *
408 * <a name="step_5-Step5setup_system"></a>
409 * <h4>Step5::setup_system</h4>
410 *
411
412 *
413 * This is the function <code>make_grid</code> from the previous
414 * example, minus the generation of the grid. Everything else is unchanged:
415 *
416 * @code
417 *   template <int dim>
418 *   void Step5<dim>::setup_system()
419 *   {
420 *   dof_handler.distribute_dofs(fe);
421 *  
422 *   std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
423 *   << std::endl;
424 *  
425 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
426 *   DoFTools::make_sparsity_pattern(dof_handler, dsp);
427 *   sparsity_pattern.copy_from(dsp);
428 *  
429 *   system_matrix.reinit(sparsity_pattern);
430 *  
431 *   solution.reinit(dof_handler.n_dofs());
432 *   system_rhs.reinit(dof_handler.n_dofs());
433 *   }
434 *  
435 *  
436 *  
437 * @endcode
438 *
439 *
440 * <a name="step_5-Step5assemble_system"></a>
441 * <h4>Step5::assemble_system</h4>
442 *
443
444 *
445 * As in the previous examples, this function is not changed much with regard
446 * to its functionality, but there are still some optimizations which we will
447 * show. For this, it is important to note that if efficient solvers are used
448 * (such as the preconditioned CG method), assembling the matrix and right hand
449 * side can take a comparable time, and you should think about using one or
450 * two optimizations at some places.
451 *
452
453 *
454 * The first parts of the function are completely unchanged from before:
455 *
456 * @code
457 *   template <int dim>
458 *   void Step5<dim>::assemble_system()
459 *   {
460 *   const QGauss<dim> quadrature_formula(fe.degree + 1);
461 *  
462 *   FEValues<dim> fe_values(fe,
463 *   quadrature_formula,
466 *  
467 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
468 *  
469 *   FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
470 *   Vector<double> cell_rhs(dofs_per_cell);
471 *  
472 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
473 *  
474 * @endcode
475 *
476 * Next is the typical loop over all cells to compute local contributions
477 * and then to transfer them into the global matrix and vector. The only
478 * change in this part, compared to @ref step_4 "step-4", is that we will use the
479 * <code>coefficient()</code> function defined above to compute the
480 * coefficient value at each quadrature point.
481 *
482 * @code
483 *   for (const auto &cell : dof_handler.active_cell_iterators())
484 *   {
485 *   fe_values.reinit(cell);
486 *  
487 *   cell_matrix = 0.;
488 *   cell_rhs = 0.;
489 *  
490 *   for (const unsigned int q_index : fe_values.quadrature_point_indices())
491 *   {
492 *   const double current_coefficient =
493 *   coefficient(fe_values.quadrature_point(q_index));
494 *   for (const unsigned int i : fe_values.dof_indices())
495 *   {
496 *   for (const unsigned int j : fe_values.dof_indices())
497 *   cell_matrix(i, j) +=
498 *   (current_coefficient * // a(x_q)
499 *   fe_values.shape_grad(i, q_index) * // grad phi_i(x_q)
500 *   fe_values.shape_grad(j, q_index) * // grad phi_j(x_q)
501 *   fe_values.JxW(q_index)); // dx
502 *  
503 *   cell_rhs(i) += (fe_values.shape_value(i, q_index) * // phi_i(x_q)
504 *   1.0 * // f(x_q)
505 *   fe_values.JxW(q_index)); // dx
506 *   }
507 *   }
508 *  
509 *  
510 *   cell->get_dof_indices(local_dof_indices);
511 *   for (const unsigned int i : fe_values.dof_indices())
512 *   {
513 *   for (const unsigned int j : fe_values.dof_indices())
514 *   system_matrix.add(local_dof_indices[i],
515 *   local_dof_indices[j],
516 *   cell_matrix(i, j));
517 *  
518 *   system_rhs(local_dof_indices[i]) += cell_rhs(i);
519 *   }
520 *   }
521 *  
522 * @endcode
523 *
524 * With the matrix so built, we use zero boundary values again:
525 *
526 * @code
527 *   std::map<types::global_dof_index, double> boundary_values;
529 *   types::boundary_id(0),
531 *   boundary_values);
532 *   MatrixTools::apply_boundary_values(boundary_values,
533 *   system_matrix,
534 *   solution,
535 *   system_rhs);
536 *   }
537 *  
538 *  
539 * @endcode
540 *
541 *
542 * <a name="step_5-Step5solve"></a>
543 * <h4>Step5::solve</h4>
544 *
545
546 *
547 * The solution process again looks mostly like in the previous
548 * examples. However, we will now use a preconditioned conjugate gradient
549 * algorithm. It is not very difficult to make this change. In fact, the only
550 * thing we have to alter is that we need an object which will act as a
551 * preconditioner. We will use SSOR (symmetric successive overrelaxation),
552 * with a relaxation factor of 1.2. For this purpose, the
553 * <code>SparseMatrix</code> class has a function which does one SSOR step,
554 * and we need to package the address of this function together with the
555 * matrix on which it should act (which is the matrix to be inverted) and the
556 * relaxation factor into one object. The <code>PreconditionSSOR</code> class
557 * does this for us. (<code>PreconditionSSOR</code> class takes a template
558 * argument denoting the matrix type it is supposed to work on. The default
559 * value is <code>SparseMatrix@<double@></code>, which is exactly what we need
560 * here, so we simply stick with the default and do not specify anything in
561 * the angle brackets.)
562 *
563
564 *
565 * Note that for the present case, SSOR doesn't really perform much better
566 * than most other preconditioners (though better than no preconditioning at
567 * all). A brief comparison of different preconditioners is presented in the
568 * Results section of the next tutorial program, @ref step_6 "step-6".
569 *
570
571 *
572 * With this, the rest of the function is trivial: instead of the
573 * <code>PreconditionIdentity</code> object we have created before, we now use
574 * the preconditioner we have declared, and the CG solver will do the rest for
575 * us:
576 *
577 * @code
578 *   template <int dim>
579 *   void Step5<dim>::solve()
580 *   {
581 *   SolverControl solver_control(1000, 1e-6 * system_rhs.l2_norm());
582 *   SolverCG<Vector<double>> solver(solver_control);
583 *  
584 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
585 *   preconditioner.initialize(system_matrix, 1.2);
586 *  
587 *   solver.solve(system_matrix, solution, system_rhs, preconditioner);
588 *  
589 *   std::cout << " " << solver_control.last_step()
590 *   << " CG iterations needed to obtain convergence." << std::endl;
591 *   }
592 *  
593 *  
594 * @endcode
595 *
596 *
597 * <a name="step_5-Step5output_resultsandsettingoutputflags"></a>
598 * <h4>Step5::output_results and setting output flags</h4>
599 *
600
601 *
602 * Writing output to a file is mostly the same as for the previous tutorial.
603 * The only difference is that we now need to construct a different filename
604 * for each refinement cycle.
605 *
606
607 *
608 * The function writes the output in VTU format, a variation of the VTK format
609 * that requires less disk space because it compresses the data. Of course,
610 * there are many other formats supported by the DataOut class if you
611 * desire to use a program for visualization that doesn't understand
612 * VTK or VTU.
613 *
614 * @code
615 *   template <int dim>
616 *   void Step5<dim>::output_results(const unsigned int cycle) const
617 *   {
618 *   DataOut<dim> data_out;
619 *  
620 *   data_out.attach_dof_handler(dof_handler);
621 *   data_out.add_data_vector(solution, "solution");
622 *  
623 *   data_out.build_patches();
624 *  
625 *   std::ofstream output("solution-" + std::to_string(cycle) + ".vtu");
626 *   data_out.write_vtu(output);
627 *   }
628 *  
629 *  
630 *  
631 * @endcode
632 *
633 *
634 * <a name="step_5-Step5run"></a>
635 * <h4>Step5::run</h4>
636 *
637
638 *
639 * The second to last thing in this program is the definition of the
640 * <code>run()</code> function. In contrast to the previous programs, we will
641 * compute on a sequence of meshes that after each iteration is globally
642 * refined. The function therefore consists of a loop over 6 cycles. In each
643 * cycle, we first print the cycle number, and then have to decide what to do
644 * with the mesh. If this is not the first cycle, we simply refine the
645 * existing mesh once globally. Before running through these cycles, however,
646 * we have to generate a mesh:
647 *
648
649 *
650 * In previous examples, we have already used some of the functions from the
651 * <code>GridGenerator</code> class. Here we would like to read a grid from a
652 * file where the cells are stored and which may originate from someone else,
653 * or may be the product of a mesh generator tool.
654 *
655
656 *
657 * In order to read a grid from a file, we generate an object of data type
658 * GridIn and associate the triangulation to it (i.e. we tell it to fill our
659 * triangulation object when we ask it to read the file). Then we open the
660 * respective file and initialize the triangulation with the data in the file :
661 *
662 * @code
663 *   template <int dim>
664 *   void Step5<dim>::run()
665 *   {
666 *   GridIn<dim> grid_in;
667 *   grid_in.attach_triangulation(triangulation);
668 *   std::ifstream input_file("circle-grid.inp");
669 * @endcode
670 *
671 * We would now like to read the file. However, the input file is only for a
672 * two-dimensional triangulation, while this function is a template for
673 * arbitrary dimension. Since this is only a demonstration program, we will
674 * not use different input files for the different dimensions, but rather
675 * quickly kill the whole program if we are not in 2d. Of course, since the
676 * main function below assumes that we are working in two dimensions we
677 * could skip this check, in this version of the program, without any ill
678 * effects.
679 *
680
681 *
682 * It turns out that perhaps 90 per cent of programming errors are invalid
683 * function parameters such as invalid array sizes, etc., so we use
684 * assertions heavily throughout deal.II to catch such mistakes. For this,
685 * the <code>Assert</code> macro is a good choice, since it makes sure that
686 * the condition which is given as first argument is valid, and if not
687 * throws an exception (its second argument) which will usually terminate
688 * the program giving information where the error occurred and what the
689 * reason was. (A longer discussion of what exactly the @p Assert macro
690 * does can be found in the @ref Exceptions "exception documentation topic".)
691 * This generally reduces the time to find programming errors
692 * dramatically and we have found assertions an invaluable means to program
693 * fast.
694 *
695
696 *
697 * On the other hand, all these checks (there are over 10,000 of them in the
698 * library at present) should not slow down the program too much if you want
699 * to do large computations. To this end, the <code>Assert</code> macro is
700 * only used in debug mode and expands to nothing if in optimized
701 * mode. Therefore, while you test your program on small problems and debug
702 * it, the assertions will tell you where the problems are. Once your
703 * program is stable, you can switch off debugging and the program will run
704 * your real computations without the assertions and at maximum speed. More
705 * precisely: turning off all the checks in the library (which prevent you
706 * from calling functions with wrong arguments, walking off of arrays, etc.)
707 * by compiling your program in optimized mode usually makes things run
708 * about four times faster. Even though optimized programs are more
709 * performant, you should always develop in debug mode since it allows
710 * the library to find lots of common programming errors automatically. For
711 * those who want to try: The way to switch from debug mode to optimized
712 * mode is to recompile your program with the command <code>make
713 * release</code>. The output of the <code>make</code> program should now
714 * indicate to you that the program is now compiled in optimized mode, and
715 * it will later also be linked to libraries that have been compiled for
716 * optimized mode. In order to switch back to debug mode, simply recompile
717 * with the command <code>make debug</code>.
718 *
719 * @code
720 *   Assert(dim == 2, ExcNotImplemented());
721 * @endcode
722 *
723 * ExcNotImplemented is a globally defined exception, which may be thrown
724 * whenever a piece of code has simply not been implemented for a case
725 * other than the condition checked in the assertion. Here, it would not
726 * be difficult to simply implement reading a *different* mesh file that
727 * contains a description of a 1d or 3d geometry, but this has not (yet)
728 * been implemented and so the exception is appropriate.
729 *
730
731 *
732 * Usually, one would like to use more specific
733 * exception classes, and particular in this case one would of course try
734 * to do something else if <code>dim</code> is not equal to two, e.g. create
735 * a grid using library functions. Aborting a program is usually not a good
736 * idea and assertions should really only be used for exceptional cases
737 * which should not occur, but might due to stupidity of the programmer,
738 * user, or someone else.
739 *
740
741 *
742 * So if we got past the assertion, we know that dim==2, and we can now
743 * actually read the grid. It is in UCD (unstructured cell data) format
744 * (though the convention is to use the suffix <code>inp</code> for UCD
745 * files):
746 *
747 * @code
748 *   grid_in.read_ucd(input_file);
749 * @endcode
750 *
751 * If you like to use another input format, you have to use one of the other
752 * <code>grid_in.read_xxx</code> function. (See the documentation of the
753 * <code>GridIn</code> class to find out what input formats are presently
754 * supported.)
755 *
756
757 *
758 * The grid in the file describes a circle. Therefore we have to use a
759 * manifold object which tells the triangulation where to put new points on
760 * the boundary when the grid is refined. Unlike @ref step_1 "step-1", since GridIn does
761 * not know that the domain has a circular boundary (unlike
762 * GridGenerator::hyper_shell) we have to explicitly attach a manifold to
763 * the boundary after creating the triangulation to get the correct result
764 * when we refine the mesh.
765 *
766 * @code
767 *   const SphericalManifold<dim> boundary;
768 *   triangulation.set_all_manifold_ids_on_boundary(0);
769 *   triangulation.set_manifold(0, boundary);
770 *  
771 *   for (unsigned int cycle = 0; cycle < 6; ++cycle)
772 *   {
773 *   std::cout << "Cycle " << cycle << ':' << std::endl;
774 *  
775 *   if (cycle != 0)
776 *   triangulation.refine_global(1);
777 *  
778 * @endcode
779 *
780 * Now that we have a mesh for sure, we write some output and do all the
781 * things that we have already seen in the previous examples.
782 *
783 * @code
784 *   std::cout << " Number of active cells: "
785 *   << triangulation.n_active_cells()
786 *   << std::endl
787 *   << " Total number of cells: "
788 *   << triangulation.n_cells()
789 *   << std::endl;
790 *  
791 *   setup_system();
792 *   assemble_system();
793 *   solve();
794 *   output_results(cycle);
795 *   }
796 *   }
797 *  
798 *  
799 * @endcode
800 *
801 *
802 * <a name="step_5-Thecodemaincodefunction"></a>
803 * <h3>The <code>main</code> function</h3>
804 *
805
806 *
807 * The main function looks mostly like the one in the previous example, so we
808 * won't comment on it further:
809 *
810 * @code
811 *   int main()
812 *   {
813 *   Step5<2> laplace_problem_2d;
814 *   laplace_problem_2d.run();
815 *   return 0;
816 *   }
817 * @endcode
818<a name="step_5-Results"></a><h1>Results</h1>
819
820
821
822Here is the console output:
823@code
824Cycle 0:
825 Number of active cells: 20
826 Total number of cells: 20
827 Number of degrees of freedom: 25
828 8 CG iterations needed to obtain convergence.
829Cycle 1:
830 Number of active cells: 80
831 Total number of cells: 100
832 Number of degrees of freedom: 89
833 12 CG iterations needed to obtain convergence.
834Cycle 2:
835 Number of active cells: 320
836 Total number of cells: 420
837 Number of degrees of freedom: 337
838 21 CG iterations needed to obtain convergence.
839Cycle 3:
840 Number of active cells: 1280
841 Total number of cells: 1700
842 Number of degrees of freedom: 1313
843 38 CG iterations needed to obtain convergence.
844Cycle 4:
845 Number of active cells: 5120
846 Total number of cells: 6820
847 Number of degrees of freedom: 5185
848 70 CG iterations needed to obtain convergence.
849Cycle 5:
850 Number of active cells: 20480
851 Total number of cells: 27300
852 Number of degrees of freedom: 20609
853 136 CG iterations needed to obtain convergence.
854@endcode
855
856
857
858In each cycle, the number of cells quadruples and the number of CG
859iterations roughly doubles.
860Also, in each cycle, the program writes one output graphic file in VTU
861format. They are depicted in the following:
862
863<table width="100%">
864 <tr>
865 <td>
866 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-0-r9.2.png" alt="">
867 </td>
868 <td>
869 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-1-r9.2.png" alt="">
870 </td>
871 </tr>
872 <tr>
873 <td>
874 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-2-r9.2.png" alt="">
875 </td>
876 <td>
877 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-3-r9.2.png" alt="">
878 </td>
879 </tr>
880 <tr>
881 <td>
882 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-4-r9.2.png" alt="">
883 </td>
884 <td>
885 <img src="https://www.dealii.org/images/steps/developer/step-5.solution-5-r9.2.png" alt="">
886 </td>
887 </tr>
888</table>
889
890
891
892Due to the variable coefficient (the curvature there is reduced by the
893same factor by which the coefficient is increased), the top region of
894the solution is flattened. The gradient of the solution is
895discontinuous along the interface, although this is not very clearly
896visible in the pictures above. We will look at this in more detail in
897the next example.
898
899The pictures also show that the solution computed by this program is
900actually pretty wrong on a very coarse mesh (its magnitude is
901wrong). That's because no numerical method guarantees that the solution
902on a coarse mesh is particularly accurate -- but we know that the
903solution <i>converges</i> to the exact solution, and indeed you can
904see how the solutions from one mesh to the next seem to not change
905very much any more at the end.
906 *
907 *
908<a name="step_5-PlainProg"></a>
909<h1> The plain program</h1>
910@include "step-5.cc"
911*/
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
Point< 2 > first
Definition grid_out.cc:4623
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())
Definition loop.h:564
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.)
Definition advection.h:74
void apply_boundary_values(const std::map< types::global_dof_index, number > &boundary_values, SparseMatrix< number > &matrix, Vector< number > &solution, Vector< number > &right_hand_side, const bool eliminate_columns=true)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation