Loading [MathJax]/extensions/TeX/newcommand.js
 Reference documentation for deal.II version 9.3.3
\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\}}
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
step-29.h
Go to the documentation of this file.
1,
408 * Vector<double> &values) const override
409 * {
410 * Assert(values.size() == 2, ExcDimensionMismatch(values.size(), 2));
411 *
412 * values(0) = 1;
413 * values(1) = 0;
414 * }
415 *
416 * virtual void
417 * vector_value_list(const std::vector<Point<dim>> &points,
418 * std::vector<Vector<double>> & value_list) const override
419 * {
420 * Assert(value_list.size() == points.size(),
421 * ExcDimensionMismatch(value_list.size(), points.size()));
422 *
423 * for (unsigned int p = 0; p < points.size(); ++p)
424 * DirichletBoundaryValues<dim>::vector_value(points[p], value_list[p]);
425 * }
426 * };
427 *
428 * @endcode
429 *
430 *
431 * <a name="ThecodeParameterReadercodeclass"></a>
432 * <h3>The <code>ParameterReader</code> class</h3>
433 *
434
435 *
436 * The next class is responsible for preparing the ParameterHandler object
437 * and reading parameters from an input file. It includes a function
438 * <code>declare_parameters</code> that declares all the necessary
439 * parameters and a <code>read_parameters</code> function that is called
440 * from outside to initiate the parameter reading process.
441 *
442 * @code
443 * class ParameterReader : public Subscriptor
444 * {
445 * public:
446 * ParameterReader(ParameterHandler &);
447 * void read_parameters(const std::string &);
448 *
449 * private:
450 * void declare_parameters();
451 * ParameterHandler &prm;
452 * };
453 *
454 * @endcode
455 *
456 * The constructor stores a reference to the ParameterHandler object that is
457 * passed to it:
458 *
459 * @code
460 * ParameterReader::ParameterReader(ParameterHandler &paramhandler)
461 * : prm(paramhandler)
462 * {}
463 *
464 * @endcode
465 *
466 *
467 * <a name="codeParameterReaderdeclare_parameterscode"></a>
468 * <h4><code>ParameterReader::declare_parameters</code></h4>
469 *
470
471 *
472 * The <code>declare_parameters</code> function declares all the parameters
473 * that our ParameterHandler object will be able to read from input files,
474 * along with their types, range conditions and the subsections they appear
475 * in. We will wrap all the entries that go into a section in a pair of
476 * braces to force the editor to indent them by one level, making it simpler
477 * to read which entries together form a section:
478 *
479 * @code
480 * void ParameterReader::declare_parameters()
481 * {
482 * @endcode
483 *
484 * Parameters for mesh and geometry include the number of global
485 * refinement steps that are applied to the initial coarse mesh and the
486 * focal distance @f$d@f$ of the transducer lens. For the number of refinement
487 * steps, we allow integer values in the range @f$[0,\infty)@f$, where the
488 * omitted second argument to the Patterns::Integer object denotes the
489 * half-open interval. For the focal distance any number greater than
490 * zero is accepted:
491 *
492 * @code
493 * prm.enter_subsection("Mesh & geometry parameters");
494 * {
495 * prm.declare_entry("Number of refinements",
496 * "6",
498 * "Number of global mesh refinement steps "
499 * "applied to initial coarse grid");
500 *
501 * prm.declare_entry("Focal distance",
502 * "0.3",
503 * Patterns::Double(0),
504 * "Distance of the focal point of the lens "
505 * "to the x-axis");
506 * }
507 * prm.leave_subsection();
508 *
509 * @endcode
510 *
511 * The next subsection is devoted to the physical parameters appearing in
512 * the equation, which are the frequency @f$\omega@f$ and wave speed
513 * @f$c@f$. Again, both need to lie in the half-open interval @f$[0,\infty)@f$
514 * represented by calling the Patterns::Double class with only the left
515 * end-point as argument:
516 *
517 * @code
518 * prm.enter_subsection("Physical constants");
519 * {
520 * prm.declare_entry("c", "1.5e5", Patterns::Double(0), "Wave speed");
521 *
522 * prm.declare_entry("omega", "5.0e7", Patterns::Double(0), "Frequency");
523 * }
524 * prm.leave_subsection();
525 *
526 *
527 * @endcode
528 *
529 * Last but not least we would like to be able to change some properties
530 * of the output, like filename and format, through entries in the
531 * configuration file, which is the purpose of the last subsection:
532 *
533 * @code
534 * prm.enter_subsection("Output parameters");
535 * {
536 * prm.declare_entry("Output filename",
537 * "solution",
539 * "Name of the output file (without extension)");
540 *
541 * @endcode
542 *
543 * Since different output formats may require different parameters for
544 * generating output (like for example, postscript output needs
545 * viewpoint angles, line widths, colors etc), it would be cumbersome if
546 * we had to declare all these parameters by hand for every possible
547 * output format supported in the library. Instead, each output format
548 * has a <code>FormatFlags::declare_parameters</code> function, which
549 * declares all the parameters specific to that format in an own
550 * subsection. The following call of
552 * <code>declare_parameters</code> for all available output formats, so
553 * that for each format an own subsection will be created with
554 * parameters declared for that particular output format. (The actual
555 * value of the template parameter in the call, <code>@<1@></code>
556 * above, does not matter here: the function does the same work
557 * independent of the dimension, but happens to be in a
558 * template-parameter-dependent class.) To find out what parameters
559 * there are for which output format, you can either consult the
560 * documentation of the DataOutBase class, or simply run this program
561 * without a parameter file present. It will then create a file with all
562 * declared parameters set to their default values, which can
563 * conveniently serve as a starting point for setting the parameters to
564 * the values you desire.
565 *
566 * @code
568 * }
569 * prm.leave_subsection();
570 * }
571 *
572 * @endcode
573 *
574 *
575 * <a name="codeParameterReaderread_parameterscode"></a>
576 * <h4><code>ParameterReader::read_parameters</code></h4>
577 *
578
579 *
580 * This is the main function in the ParameterReader class. It gets called
581 * from outside, first declares all the parameters, and then reads them from
582 * the input file whose filename is provided by the caller. After the call
583 * to this function is complete, the <code>prm</code> object can be used to
584 * retrieve the values of the parameters read in from the file :
585 *
586 * @code
587 * void ParameterReader::read_parameters(const std::string &parameter_file)
588 * {
589 * declare_parameters();
590 *
591 * prm.parse_input(parameter_file);
592 * }
593 *
594 *
595 *
596 * @endcode
597 *
598 *
599 * <a name="ThecodeComputeIntensitycodeclass"></a>
600 * <h3>The <code>ComputeIntensity</code> class</h3>
601 *
602
603 *
604 * As mentioned in the introduction, the quantity that we are really after
605 * is the spatial distribution of the intensity of the ultrasound wave,
606 * which corresponds to @f$|u|=\sqrt{v^2+w^2}@f$. Now we could just be content
607 * with having @f$v@f$ and @f$w@f$ in our output, and use a suitable visualization
608 * or postprocessing tool to derive @f$|u|@f$ from the solution we
609 * computed. However, there is also a way to output data derived from the
610 * solution in deal.II, and we are going to make use of this mechanism here.
611 *
612
613 *
614 * So far we have always used the DataOut::add_data_vector function to add
615 * vectors containing output data to a DataOut object. There is a special
616 * version of this function that in addition to the data vector has an
617 * additional argument of type DataPostprocessor. What happens when this
618 * function is used for output is that at each point where output data is to
619 * be generated, the DataPostprocessor::evaluate_scalar_field() or
620 * DataPostprocessor::evaluate_vector_field() function of the
621 * specified DataPostprocessor object is invoked to compute the output
622 * quantities from the values, the gradients and the second derivatives of
623 * the finite element function represented by the data vector (in the case
624 * of face related data, normal vectors are available as well). Hence, this
625 * allows us to output any quantity that can locally be derived from the
626 * values of the solution and its derivatives. Of course, the ultrasound
627 * intensity @f$|u|@f$ is such a quantity and its computation doesn't even
628 * involve any derivatives of @f$v@f$ or @f$w@f$.
629 *
630
631 *
632 * In practice, the DataPostprocessor class only provides an interface to
633 * this functionality, and we need to derive our own class from it in order
634 * to implement the functions specified by the interface. In the most
635 * general case one has to implement several member functions but if the
636 * output quantity is a single scalar then some of this boilerplate code can
637 * be handled by a more specialized class, DataPostprocessorScalar and we
638 * can derive from that one instead. This is what the
639 * <code>ComputeIntensity</code> class does:
640 *
641 * @code
642 * template <int dim>
643 * class ComputeIntensity : public DataPostprocessorScalar<dim>
644 * {
645 * public:
646 * ComputeIntensity();
647 *
648 * virtual void evaluate_vector_field(
650 * std::vector<Vector<double>> &computed_quantities) const override;
651 * };
652 *
653 * @endcode
654 *
655 * In the constructor, we need to call the constructor of the base class
656 * with two arguments. The first denotes the name by which the single scalar
657 * quantity computed by this class should be represented in output files. In
658 * our case, the postprocessor has @f$|u|@f$ as output, so we use "Intensity".
659 *
660
661 *
662 * The second argument is a set of flags that indicate which data is needed
663 * by the postprocessor in order to compute the output quantities. This can
664 * be any subset of update_values, update_gradients and update_hessians
665 * (and, in the case of face data, also update_normal_vectors), which are
666 * documented in UpdateFlags. Of course, computation of the derivatives
667 * requires additional resources, so only the flags for data that are really
668 * needed should be given here, just as we do when we use FEValues objects.
669 * In our case, only the function values of @f$v@f$ and @f$w@f$ are needed to
670 * compute @f$|u|@f$, so we're good with the update_values flag.
671 *
672 * @code
673 * template <int dim>
674 * ComputeIntensity<dim>::ComputeIntensity()
675 * : DataPostprocessorScalar<dim>("Intensity", update_values)
676 * {}
677 *
678 *
679 * @endcode
680 *
681 * The actual postprocessing happens in the following function. Its input is
682 * an object that stores values of the function (which is here vector-valued)
683 * representing the data vector given to DataOut::add_data_vector, evaluated
684 * at all evaluation points where we generate output, and some tensor objects
685 * representing derivatives (that we don't use here since @f$|u|@f$ is computed
686 * from just @f$v@f$ and @f$w@f$). The derived quantities are returned in the
687 * <code>computed_quantities</code> vector. Remember that this function may
688 * only use data for which the respective update flag is specified by
689 * <code>get_needed_update_flags</code>. For example, we may not use the
690 * derivatives here, since our implementation of
691 * <code>get_needed_update_flags</code> requests that only function values
692 * are provided.
693 *
694 * @code
695 * template <int dim>
696 * void ComputeIntensity<dim>::evaluate_vector_field(
698 * std::vector<Vector<double>> & computed_quantities) const
699 * {
700 * Assert(computed_quantities.size() == inputs.solution_values.size(),
701 * ExcDimensionMismatch(computed_quantities.size(),
702 * inputs.solution_values.size()));
703 *
704 * @endcode
705 *
706 * The computation itself is straightforward: We iterate over each
707 * entry in the output vector and compute @f$|u|@f$ from the
708 * corresponding values of @f$v@f$ and @f$w@f$. We do this by creating a
709 * complex number @f$u@f$ and then calling `std::abs()` on the
710 * result. (One may be tempted to call `std::norm()`, but in a
711 * historical quirk, the C++ committee decided that `std::norm()`
712 * should return the <i>square</i> of the absolute value --
713 * thereby not satisfying the properties mathematicians require of
714 * something called a "norm".)
715 *
716 * @code
717 * for (unsigned int i = 0; i < computed_quantities.size(); i++)
718 * {
719 * Assert(computed_quantities[i].size() == 1,
720 * ExcDimensionMismatch(computed_quantities[i].size(), 1));
721 * Assert(inputs.solution_values[i].size() == 2,
722 * ExcDimensionMismatch(inputs.solution_values[i].size(), 2));
723 *
724 * const std::complex<double> u(inputs.solution_values[i](0),
725 * inputs.solution_values[i](1));
726 *
727 * computed_quantities[i](0) = std::abs(u);
728 * }
729 * }
730 *
731 *
732 * @endcode
733 *
734 *
735 * <a name="ThecodeUltrasoundProblemcodeclass"></a>
736 * <h3>The <code>UltrasoundProblem</code> class</h3>
737 *
738
739 *
740 * Finally here is the main class of this program. It's member functions
741 * are very similar to the previous examples, in particular @ref step_4 "step-4", and the
742 * list of member variables does not contain any major surprises either.
743 * The ParameterHandler object that is passed to the constructor is stored
744 * as a reference to allow easy access to the parameters from all functions
745 * of the class. Since we are working with vector valued finite elements,
746 * the FE object we are using is of type FESystem.
747 *
748 * @code
749 * template <int dim>
750 * class UltrasoundProblem
751 * {
752 * public:
753 * UltrasoundProblem(ParameterHandler &);
754 * void run();
755 *
756 * private:
757 * void make_grid();
758 * void setup_system();
759 * void assemble_system();
760 * void solve();
761 * void output_results() const;
762 *
763 * ParameterHandler &prm;
764 *
766 * DoFHandler<dim> dof_handler;
767 * FESystem<dim> fe;
768 *
769 * SparsityPattern sparsity_pattern;
770 * SparseMatrix<double> system_matrix;
771 * Vector<double> solution, system_rhs;
772 * };
773 *
774 *
775 *
776 * @endcode
777 *
778 * The constructor takes the ParameterHandler object and stores it in a
779 * reference. It also initializes the DoF-Handler and the finite element
780 * system, which consists of two copies of the scalar Q1 field, one for @f$v@f$
781 * and one for @f$w@f$:
782 *
783 * @code
784 * template <int dim>
785 * UltrasoundProblem<dim>::UltrasoundProblem(ParameterHandler &param)
786 * : prm(param)
787 * , dof_handler(triangulation)
788 * , fe(FE_Q<dim>(1), 2)
789 * {}
790 *
791 * @endcode
792 *
793 *
794 * <a name="codeUltrasoundProblemmake_gridcode"></a>
795 * <h4><code>UltrasoundProblem::make_grid</code></h4>
796 *
797
798 *
799 * Here we setup the grid for our domain. As mentioned in the exposition,
800 * the geometry is just a unit square (in 2d) with the part of the boundary
801 * that represents the transducer lens replaced by a sector of a circle.
802 *
803 * @code
804 * template <int dim>
805 * void UltrasoundProblem<dim>::make_grid()
806 * {
807 * @endcode
808 *
809 * First we generate some logging output and start a timer so we can
810 * compute execution time when this function is done:
811 *
812 * @code
813 * std::cout << "Generating grid... ";
814 * Timer timer;
815 *
816 * @endcode
817 *
818 * Then we query the values for the focal distance of the transducer lens
819 * and the number of mesh refinement steps from our ParameterHandler
820 * object:
821 *
822 * @code
823 * prm.enter_subsection("Mesh & geometry parameters");
824 *
825 * const double focal_distance = prm.get_double("Focal distance");
826 * const unsigned int n_refinements = prm.get_integer("Number of refinements");
827 *
828 * prm.leave_subsection();
829 *
830 * @endcode
831 *
832 * Next, two points are defined for position and focal point of the
833 * transducer lens, which is the center of the circle whose segment will
834 * form the transducer part of the boundary. Notice that this is the only
835 * point in the program where things are slightly different in 2D and 3D.
836 * Even though this tutorial only deals with the 2D case, the necessary
837 * additions to make this program functional in 3D are so minimal that we
838 * opt for including them:
839 *
840 * @code
841 * const Point<dim> transducer =
842 * (dim == 2) ? Point<dim>(0.5, 0.0) : Point<dim>(0.5, 0.5, 0.0);
843 * const Point<dim> focal_point = (dim == 2) ?
844 * Point<dim>(0.5, focal_distance) :
845 * Point<dim>(0.5, 0.5, focal_distance);
846 *
847 *
848 * @endcode
849 *
850 * As initial coarse grid we take a simple unit square with 5 subdivisions
851 * in each direction. The number of subdivisions is chosen so that the
852 * line segment @f$[0.4,0.6]@f$ that we want to designate as the transducer
853 * boundary is spanned by a single face. Then we step through all cells to
854 * find the faces where the transducer is to be located, which in fact is
855 * just the single edge from 0.4 to 0.6 on the x-axis. This is where we
856 * want the refinements to be made according to a circle shaped boundary,
857 * so we mark this edge with a different manifold indicator. Since we will
858 * Dirichlet boundary conditions on the transducer, we also change its
859 * boundary indicator.
860 *
861 * @code
863 *
864 * for (auto &cell : triangulation.cell_iterators())
865 * for (const auto &face : cell->face_iterators())
866 * if (face->at_boundary() &&
867 * ((face->center() - transducer).norm_square() < 0.01))
868 * {
869 * face->set_boundary_id(1);
870 * face->set_manifold_id(1);
871 * }
872 * @endcode
873 *
874 * For the circle part of the transducer lens, a SphericalManifold object
875 * is used (which, of course, in 2D just represents a circle), with center
876 * computed as above.
877 *
878 * @code
879 * triangulation.set_manifold(1, SphericalManifold<dim>(focal_point));
880 *
881 * @endcode
882 *
883 * Now global refinement is executed. Cells near the transducer location
884 * will be automatically refined according to the circle shaped boundary
885 * of the transducer lens:
886 *
887 * @code
888 * triangulation.refine_global(n_refinements);
889 *
890 * @endcode
891 *
892 * Lastly, we generate some more logging output. We stop the timer and
893 * query the number of CPU seconds elapsed since the beginning of the
894 * function:
895 *
896 * @code
897 * timer.stop();
898 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
899 *
900 * std::cout << " Number of active cells: " << triangulation.n_active_cells()
901 * << std::endl;
902 * }
903 *
904 *
905 * @endcode
906 *
907 *
908 * <a name="codeUltrasoundProblemsetup_systemcode"></a>
909 * <h4><code>UltrasoundProblem::setup_system</code></h4>
910 *
911
912 *
913 * Initialization of the system matrix, sparsity patterns and vectors are
914 * the same as in previous examples and therefore do not need further
915 * comment. As in the previous function, we also output the run time of what
916 * we do here:
917 *
918 * @code
919 * template <int dim>
920 * void UltrasoundProblem<dim>::setup_system()
921 * {
922 * std::cout << "Setting up system... ";
923 * Timer timer;
924 *
925 * dof_handler.distribute_dofs(fe);
926 *
927 * DynamicSparsityPattern dsp(dof_handler.n_dofs(), dof_handler.n_dofs());
928 * DoFTools::make_sparsity_pattern(dof_handler, dsp);
929 * sparsity_pattern.copy_from(dsp);
930 *
931 * system_matrix.reinit(sparsity_pattern);
932 * system_rhs.reinit(dof_handler.n_dofs());
933 * solution.reinit(dof_handler.n_dofs());
934 *
935 * timer.stop();
936 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
937 *
938 * std::cout << " Number of degrees of freedom: " << dof_handler.n_dofs()
939 * << std::endl;
940 * }
941 *
942 *
943 * @endcode
944 *
945 *
946 * <a name="codeUltrasoundProblemassemble_systemcode"></a>
947 * <h4><code>UltrasoundProblem::assemble_system</code></h4>
948 *
949
950 *
951 * As before, this function takes care of assembling the system matrix and
952 * right hand side vector:
953 *
954 * @code
955 * template <int dim>
956 * void UltrasoundProblem<dim>::assemble_system()
957 * {
958 * std::cout << "Assembling system matrix... ";
959 * Timer timer;
960 *
961 * @endcode
962 *
963 * First we query wavespeed and frequency from the ParameterHandler object
964 * and store them in local variables, as they will be used frequently
965 * throughout this function.
966 *
967
968 *
969 *
970 * @code
971 * prm.enter_subsection("Physical constants");
972 *
973 * const double omega = prm.get_double("omega"), c = prm.get_double("c");
974 *
975 * prm.leave_subsection();
976 *
977 * @endcode
978 *
979 * As usual, for computing integrals ordinary Gauss quadrature rule is
980 * used. Since our bilinear form involves boundary integrals on
981 * @f$\Gamma_2@f$, we also need a quadrature rule for surface integration on
982 * the faces, which are @f$dim-1@f$ dimensional:
983 *
984 * @code
985 * QGauss<dim> quadrature_formula(fe.degree + 1);
986 * QGauss<dim - 1> face_quadrature_formula(fe.degree + 1);
987 *
988 * const unsigned int n_q_points = quadrature_formula.size(),
989 * n_face_q_points = face_quadrature_formula.size(),
990 * dofs_per_cell = fe.n_dofs_per_cell();
991 *
992 * @endcode
993 *
994 * The FEValues objects will evaluate the shape functions for us. For the
995 * part of the bilinear form that involves integration on @f$\Omega@f$, we'll
996 * need the values and gradients of the shape functions, and of course the
997 * quadrature weights. For the terms involving the boundary integrals,
998 * only shape function values and the quadrature weights are necessary.
999 *
1000 * @code
1001 * FEValues<dim> fe_values(fe,
1002 * quadrature_formula,
1003 * update_values | update_gradients |
1004 * update_JxW_values);
1005 *
1006 * FEFaceValues<dim> fe_face_values(fe,
1007 * face_quadrature_formula,
1008 * update_values | update_JxW_values);
1009 *
1010 * @endcode
1011 *
1012 * As usual, the system matrix is assembled cell by cell, and we need a
1013 * matrix for storing the local cell contributions as well as an index
1014 * vector to transfer the cell contributions to the appropriate location
1015 * in the global system matrix after.
1016 *
1017 * @code
1018 * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
1019 * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1020 *
1021 * for (const auto &cell : dof_handler.active_cell_iterators())
1022 * {
1023 * @endcode
1024 *
1025 * On each cell, we first need to reset the local contribution matrix
1026 * and request the FEValues object to compute the shape functions for
1027 * the current cell:
1028 *
1029 * @code
1030 * cell_matrix = 0;
1031 * fe_values.reinit(cell);
1032 *
1033 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1034 * {
1035 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1036 * {
1037 * @endcode
1038 *
1039 * At this point, it is important to keep in mind that we are
1040 * dealing with a finite element system with two
1041 * components. Due to the way we constructed this FESystem,
1042 * namely as the Cartesian product of two scalar finite
1043 * element fields, each shape function has only a single
1044 * nonzero component (they are, in deal.II lingo, @ref
1045 * GlossPrimitive "primitive"). Hence, each shape function
1046 * can be viewed as one of the @f$\phi@f$'s or @f$\psi@f$'s from the
1047 * introduction, and similarly the corresponding degrees of
1048 * freedom can be attributed to either @f$\alpha@f$ or @f$\beta@f$.
1049 * As we iterate through all the degrees of freedom on the
1050 * current cell however, they do not come in any particular
1051 * order, and so we cannot decide right away whether the DoFs
1052 * with index @f$i@f$ and @f$j@f$ belong to the real or imaginary part
1053 * of our solution. On the other hand, if you look at the
1054 * form of the system matrix in the introduction, this
1055 * distinction is crucial since it will determine to which
1056 * block in the system matrix the contribution of the current
1057 * pair of DoFs will go and hence which quantity we need to
1058 * compute from the given two shape functions. Fortunately,
1059 * the FESystem object can provide us with this information,
1060 * namely it has a function
1061 * FESystem::system_to_component_index, that for each local
1062 * DoF index returns a pair of integers of which the first
1063 * indicates to which component of the system the DoF
1064 * belongs. The second integer of the pair indicates which
1065 * index the DoF has in the scalar base finite element field,
1066 * but this information is not relevant here. If you want to
1067 * know more about this function and the underlying scheme
1068 * behind primitive vector valued elements, take a look at
1069 * @ref step_8 "step-8" or the @ref vector_valued module, where these topics
1070 * are explained in depth.
1071 *
1072 * @code
1073 * if (fe.system_to_component_index(i).first ==
1074 * fe.system_to_component_index(j).first)
1075 * {
1076 * @endcode
1077 *
1078 * If both DoFs @f$i@f$ and @f$j@f$ belong to same component,
1079 * i.e. their shape functions are both @f$\phi@f$'s or both
1080 * @f$\psi@f$'s, the contribution will end up in one of the
1081 * diagonal blocks in our system matrix, and since the
1082 * corresponding entries are computed by the same formula,
1083 * we do not bother if they actually are @f$\phi@f$ or @f$\psi@f$
1084 * shape functions. We can simply compute the entry by
1085 * iterating over all quadrature points and adding up
1086 * their contributions, where values and gradients of the
1087 * shape functions are supplied by our FEValues object.
1088 *
1089
1090 *
1091 *
1092 * @code
1093 * for (unsigned int q_point = 0; q_point < n_q_points;
1094 * ++q_point)
1095 * cell_matrix(i, j) +=
1096 * (((fe_values.shape_value(i, q_point) *
1097 * fe_values.shape_value(j, q_point)) *
1098 * (-omega * omega) +
1099 * (fe_values.shape_grad(i, q_point) *
1100 * fe_values.shape_grad(j, q_point)) *
1101 * c * c) *
1102 * fe_values.JxW(q_point));
1103 *
1104 * @endcode
1105 *
1106 * You might think that we would have to specify which
1107 * component of the shape function we'd like to evaluate
1108 * when requesting shape function values or gradients from
1109 * the FEValues object. However, as the shape functions
1110 * are primitive, they have only one nonzero component,
1111 * and the FEValues class is smart enough to figure out
1112 * that we are definitely interested in this one nonzero
1113 * component.
1114 *
1115 * @code
1116 * }
1117 * }
1118 * }
1119 *
1120 *
1121 * @endcode
1122 *
1123 * We also have to add contributions due to boundary terms. To this
1124 * end, we loop over all faces of the current cell and see if first it
1125 * is at the boundary, and second has the correct boundary indicator
1126 * associated with @f$\Gamma_2@f$, the part of the boundary where we have
1127 * absorbing boundary conditions:
1128 *
1129 * @code
1130 * for (const auto face_no : cell->face_indices())
1131 * if (cell->face(face_no)->at_boundary() &&
1132 * (cell->face(face_no)->boundary_id() == 0))
1133 * {
1134 * @endcode
1135 *
1136 * These faces will certainly contribute to the off-diagonal
1137 * blocks of the system matrix, so we ask the FEFaceValues
1138 * object to provide us with the shape function values on this
1139 * face:
1140 *
1141 * @code
1142 * fe_face_values.reinit(cell, face_no);
1143 *
1144 *
1145 * @endcode
1146 *
1147 * Next, we loop through all DoFs of the current cell to find
1148 * pairs that belong to different components and both have
1149 * support on the current face_no:
1150 *
1151 * @code
1152 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1153 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1154 * if ((fe.system_to_component_index(i).first !=
1155 * fe.system_to_component_index(j).first) &&
1156 * fe.has_support_on_face(i, face_no) &&
1157 * fe.has_support_on_face(j, face_no))
1158 * @endcode
1159 *
1160 * The check whether shape functions have support on a
1161 * face is not strictly necessary: if we don't check for
1162 * it we would simply add up terms to the local cell
1163 * matrix that happen to be zero because at least one of
1164 * the shape functions happens to be zero. However, we can
1165 * save that work by adding the checks above.
1166 *
1167
1168 *
1169 * In either case, these DoFs will contribute to the
1170 * boundary integrals in the off-diagonal blocks of the
1171 * system matrix. To compute the integral, we loop over
1172 * all the quadrature points on the face and sum up the
1173 * contribution weighted with the quadrature weights that
1174 * the face quadrature rule provides. In contrast to the
1175 * entries on the diagonal blocks, here it does matter
1176 * which one of the shape functions is a @f$\psi@f$ and which
1177 * one is a @f$\phi@f$, since that will determine the sign of
1178 * the entry. We account for this by a simple conditional
1179 * statement that determines the correct sign. Since we
1180 * already checked that DoF @f$i@f$ and @f$j@f$ belong to
1181 * different components, it suffices here to test for one
1182 * of them to which component it belongs.
1183 *
1184 * @code
1185 * for (unsigned int q_point = 0; q_point < n_face_q_points;
1186 * ++q_point)
1187 * cell_matrix(i, j) +=
1188 * ((fe.system_to_component_index(i).first == 0) ? -1 :
1189 * 1) *
1190 * fe_face_values.shape_value(i, q_point) *
1191 * fe_face_values.shape_value(j, q_point) * c * omega *
1192 * fe_face_values.JxW(q_point);
1193 * }
1194 *
1195 * @endcode
1196 *
1197 * Now we are done with this cell and have to transfer its
1198 * contributions from the local to the global system matrix. To this
1199 * end, we first get a list of the global indices of the this cells
1200 * DoFs...
1201 *
1202 * @code
1203 * cell->get_dof_indices(local_dof_indices);
1204 *
1205 *
1206 * @endcode
1207 *
1208 * ...and then add the entries to the system matrix one by one:
1209 *
1210 * @code
1211 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
1212 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
1213 * system_matrix.add(local_dof_indices[i],
1214 * local_dof_indices[j],
1215 * cell_matrix(i, j));
1216 * }
1217 *
1218 *
1219 * @endcode
1220 *
1221 * The only thing left are the Dirichlet boundary values on @f$\Gamma_1@f$,
1222 * which is characterized by the boundary indicator 1. The Dirichlet
1223 * values are provided by the <code>DirichletBoundaryValues</code> class
1224 * we defined above:
1225 *
1226 * @code
1227 * std::map<types::global_dof_index, double> boundary_values;
1228 * VectorTools::interpolate_boundary_values(dof_handler,
1229 * 1,
1230 * DirichletBoundaryValues<dim>(),
1231 * boundary_values);
1232 *
1233 * MatrixTools::apply_boundary_values(boundary_values,
1234 * system_matrix,
1235 * solution,
1236 * system_rhs);
1237 *
1238 * timer.stop();
1239 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1240 * }
1241 *
1242 *
1243 *
1244 * @endcode
1245 *
1246 *
1247 * <a name="codeUltrasoundProblemsolvecode"></a>
1248 * <h4><code>UltrasoundProblem::solve</code></h4>
1249 *
1250
1251 *
1252 * As already mentioned in the introduction, the system matrix is neither
1253 * symmetric nor definite, and so it is not quite obvious how to come up
1254 * with an iterative solver and a preconditioner that do a good job on this
1255 * matrix. We chose instead to go a different way and solve the linear
1256 * system with the sparse LU decomposition provided by UMFPACK. This is
1257 * often a good first choice for 2D problems and works reasonably well even
1258 * for a large number of DoFs. The deal.II interface to UMFPACK is given by
1259 * the SparseDirectUMFPACK class, which is very easy to use and allows us to
1260 * solve our linear system with just 3 lines of code.
1261 *
1262
1263 *
1264 * Note again that for compiling this example program, you need to have the
1265 * deal.II library built with UMFPACK support.
1266 *
1267 * @code
1268 * template <int dim>
1269 * void UltrasoundProblem<dim>::solve()
1270 * {
1271 * std::cout << "Solving linear system... ";
1272 * Timer timer;
1273 *
1274 * @endcode
1275 *
1276 * The code to solve the linear system is short: First, we allocate an
1277 * object of the right type. The following <code>initialize</code> call
1278 * provides the matrix that we would like to invert to the
1279 * SparseDirectUMFPACK object, and at the same time kicks off the
1280 * LU-decomposition. Hence, this is also the point where most of the
1281 * computational work in this program happens.
1282 *
1283 * @code
1284 * SparseDirectUMFPACK A_direct;
1285 * A_direct.initialize(system_matrix);
1286 *
1287 * @endcode
1288 *
1289 * After the decomposition, we can use <code>A_direct</code> like a matrix
1290 * representing the inverse of our system matrix, so to compute the
1291 * solution we just have to multiply with the right hand side vector:
1292 *
1293 * @code
1294 * A_direct.vmult(solution, system_rhs);
1295 *
1296 * timer.stop();
1297 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1298 * }
1299 *
1300 *
1301 *
1302 * @endcode
1303 *
1304 *
1305 * <a name="codeUltrasoundProblemoutput_resultscode"></a>
1306 * <h4><code>UltrasoundProblem::output_results</code></h4>
1307 *
1308
1309 *
1310 * Here we output our solution @f$v@f$ and @f$w@f$ as well as the derived quantity
1311 * @f$|u|@f$ in the format specified in the parameter file. Most of the work for
1312 * deriving @f$|u|@f$ from @f$v@f$ and @f$w@f$ was already done in the implementation of
1313 * the <code>ComputeIntensity</code> class, so that the output routine is
1314 * rather straightforward and very similar to what is done in the previous
1315 * tutorials.
1316 *
1317 * @code
1318 * template <int dim>
1319 * void UltrasoundProblem<dim>::output_results() const
1320 * {
1321 * std::cout << "Generating output... ";
1322 * Timer timer;
1323 *
1324 * @endcode
1325 *
1326 * Define objects of our <code>ComputeIntensity</code> class and a DataOut
1327 * object:
1328 *
1329 * @code
1330 * ComputeIntensity<dim> intensities;
1331 * DataOut<dim> data_out;
1332 *
1333 * data_out.attach_dof_handler(dof_handler);
1334 *
1335 * @endcode
1336 *
1337 * Next we query the output-related parameters from the ParameterHandler.
1338 * The DataOut::parse_parameters call acts as a counterpart to the
1339 * DataOutInterface<1>::declare_parameters call in
1340 * <code>ParameterReader::declare_parameters</code>. It collects all the
1341 * output format related parameters from the ParameterHandler and sets the
1342 * corresponding properties of the DataOut object accordingly.
1343 *
1344 * @code
1345 * prm.enter_subsection("Output parameters");
1346 *
1347 * const std::string output_filename = prm.get("Output filename");
1348 * data_out.parse_parameters(prm);
1349 *
1350 * prm.leave_subsection();
1351 *
1352 * @endcode
1353 *
1354 * Now we put together the filename from the base name provided by the
1355 * ParameterHandler and the suffix which is provided by the DataOut class
1356 * (the default suffix is set to the right type that matches the one set
1357 * in the .prm file through parse_parameters()):
1358 *
1359 * @code
1360 * const std::string filename = output_filename + data_out.default_suffix();
1361 *
1362 * std::ofstream output(filename);
1363 *
1364 * @endcode
1365 *
1366 * The solution vectors @f$v@f$ and @f$w@f$ are added to the DataOut object in the
1367 * usual way:
1368 *
1369 * @code
1370 * std::vector<std::string> solution_names;
1371 * solution_names.emplace_back("Re_u");
1372 * solution_names.emplace_back("Im_u");
1373 *
1374 * data_out.add_data_vector(solution, solution_names);
1375 *
1376 * @endcode
1377 *
1378 * For the intensity, we just call <code>add_data_vector</code> again, but
1379 * this with our <code>ComputeIntensity</code> object as the second
1380 * argument, which effectively adds @f$|u|@f$ to the output data:
1381 *
1382 * @code
1383 * data_out.add_data_vector(solution, intensities);
1384 *
1385 * @endcode
1386 *
1387 * The last steps are as before. Note that the actual output format is now
1388 * determined by what is stated in the input file, i.e. one can change the
1389 * output format without having to re-compile this program:
1390 *
1391 * @code
1392 * data_out.build_patches();
1393 * data_out.write(output);
1394 *
1395 * timer.stop();
1396 * std::cout << "done (" << timer.cpu_time() << "s)" << std::endl;
1397 * }
1398 *
1399 *
1400 *
1401 * @endcode
1402 *
1403 *
1404 * <a name="codeUltrasoundProblemruncode"></a>
1405 * <h4><code>UltrasoundProblem::run</code></h4>
1406 *
1407
1408 *
1409 * Here we simply execute our functions one after the other:
1410 *
1411 * @code
1412 * template <int dim>
1413 * void UltrasoundProblem<dim>::run()
1414 * {
1415 * make_grid();
1416 * setup_system();
1417 * assemble_system();
1418 * solve();
1419 * output_results();
1420 * }
1421 * } // namespace Step29
1422 *
1423 *
1424 * @endcode
1425 *
1426 *
1427 * <a name="Thecodemaincodefunction"></a>
1428 * <h4>The <code>main</code> function</h4>
1429 *
1430
1431 *
1432 * Finally the <code>main</code> function of the program. It has the same
1433 * structure as in almost all of the other tutorial programs. The only
1434 * exception is that we define ParameterHandler and
1435 * <code>ParameterReader</code> objects, and let the latter read in the
1436 * parameter values from a textfile called <code>@ref step_29 "step-29".prm</code>. The
1437 * values so read are then handed over to an instance of the UltrasoundProblem
1438 * class:
1439 *
1440 * @code
1441 * int main()
1442 * {
1443 * try
1444 * {
1445 * using namespace dealii;
1446 * using namespace Step29;
1447 *
1448 * ParameterHandler prm;
1449 * ParameterReader param(prm);
1450 * param.read_parameters("step-29.prm");
1451 *
1452 * UltrasoundProblem<2> ultrasound_problem(prm);
1453 * ultrasound_problem.run();
1454 * }
1455 * catch (std::exception &exc)
1456 * {
1457 * std::cerr << std::endl
1458 * << std::endl
1459 * << "----------------------------------------------------"
1460 * << std::endl;
1461 * std::cerr << "Exception on processing: " << std::endl
1462 * << exc.what() << std::endl
1463 * << "Aborting!" << std::endl
1464 * << "----------------------------------------------------"
1465 * << std::endl;
1466 * return 1;
1467 * }
1468 * catch (...)
1469 * {
1470 * std::cerr << std::endl
1471 * << std::endl
1472 * << "----------------------------------------------------"
1473 * << std::endl;
1474 * std::cerr << "Unknown exception!" << std::endl
1475 * << "Aborting!" << std::endl
1476 * << "----------------------------------------------------"
1477 * << std::endl;
1478 * return 1;
1479 * }
1480 * return 0;
1481 * }
1482 * @endcode
1483<a name="Results"></a>
1484<a name="Results"></a><h1>Results</h1>
1485
1486
1487The current program reads its run-time parameters from an input file
1488called <code>step-29.prm</code> that looks like this:
1489@code
1490subsection Mesh & geometry parameters
1491 # Distance of the focal point of the lens to the x-axis
1492 set Focal distance = 0.3
1493
1494 # Number of global mesh refinement steps applied to initial coarse grid
1495 set Number of refinements = 5
1496end
1497
1498
1499subsection Physical constants
1500 # Wave speed
1501 set c = 1.5e5
1502
1503 # Frequency
1504 set omega = 3.0e7
1505end
1506
1507
1508subsection Output parameters
1509 # Name of the output file (without extension)
1510 set Output file = solution
1511
1512 # A name for the output format to be used
1513 set Output format = vtu
1514end
1515@endcode
1516
1517As can be seen, we set
1518@f$d=0.3@f$, which amounts to a focus of the transducer lens
1519at @f$x=0.5@f$, @f$y=0.3@f$. The coarse mesh is refined 5 times,
1520resulting in 160x160 cells, and the output is written in vtu
1521format. The parameter reader understands many more parameters
1522pertaining in particular to the generation of output, but we
1523need none of these parameters here and therefore stick with
1524their default values.
1525
1526Here's the console output of the program in debug mode:
1527
1528@code
1529> make run
1530[ 66%] Built target @ref step_29 "step-29"
1531[100%] Run @ref step_29 "step-29" with Debug configuration
1532Generating grid... done (0.820449s)
1533 Number of active cells: 25600
1534Setting up system... done (1.18392s)
1535 Number of degrees of freedom: 51842
1536Assembling system matrix... done (2.33291s)
1537Solving linear system... done (1.34837s)
1538Generating output... done (2.05782s)
1539[100%] Built target run
1540@endcode
1541
1542(Of course, execution times will differ if you run the program
1543locally.) The fact that most of the time is spent on assembling
1544the system matrix and generating output is due to the many assertions
1545that need to be checked in debug mode. In release mode these parts
1546of the program run much faster whereas solving the linear system is
1547hardly sped up at all:
1548
1549@code
1550> make run
1551[ 66%] Built target @ref step_29 "step-29"
1552Scanning dependencies of target run
1553[100%] Run @ref step_29 "step-29" with Release configuration
1554DEAL::Generating grid... done (0.0144960s)
1555DEAL:: Number of active cells: 25600
1556DEAL::Setting up system... done (0.0356880s)
1557DEAL:: Number of degrees of freedom: 51842
1558DEAL::Assembling system matrix... done (0.0436570s)
1559DEAL::Solving linear system... done (1.54733s)
1560DEAL::Generating output... done (0.720528s)
1561[100%] Built target run
1562@endcode
1563
1564The graphical output of the program looks as follows:
1565
1566
1567<table align="center" class="doxtable">
1568 <tr>
1569 <td>
1570 <img src="https://www.dealii.org/images/steps/developer/step-29.v.png" alt="v = Re(u)">
1571 </td>
1572 <td>
1573 <img src="https://www.dealii.org/images/steps/developer/step-29.w.png" alt="w = Im(u)">
1574 </td>
1575 </tr>
1576 <tr>
1577 <td colspan="2">
1578 <img src="https://www.dealii.org/images/steps/developer/step-29.intensity.png" alt="|u|">
1579 </td>
1580 </tr>
1581</table>
1582
1583The first two pictures show the real and imaginary parts of
1584@f$u@f$, whereas the last shows the intensity @f$|u|@f$. One can clearly
1585see that the intensity is focused around the focal point of the
1586lens (0.5, 0.3), and that the focus
1587is rather sharp in @f$x@f$-direction but more blurred in @f$y@f$-direction, which is a
1588consequence of the geometry of the focusing lens, its finite aperture,
1589and the wave nature of the problem.
1590
1591Because colorful graphics are always fun, and to stress the focusing
1592effects some more, here is another set of images highlighting how well
1593the intensity is actually focused in @f$x@f$-direction:
1594
1595<table align="center" class="doxtable">
1596 <tr>
1597 <td>
1598 <img src="https://www.dealii.org/images/steps/developer/step-29.surface.png" alt="|u|">
1599 </td>
1600 <td>
1601 <img src="https://www.dealii.org/images/steps/developer/step-29.contours.png" alt="|u|">
1602 </td>
1603 </tr>
1604</table>
1605
1606
1607As a final note, the structure of the program makes it easy to
1608determine which parts of the program scale nicely as the mesh is
1609refined and which parts don't. Here are the run times for 5, 6, and 7
1610global refinements:
1611
1612@code
1613> make run
1614[ 66%] Built target @ref step_29 "step-29"
1615[100%] Run @ref step_29 "step-29" with Release configuration
1616DEAL::Generating grid... done (0.0135260s)
1617DEAL:: Number of active cells: 25600
1618DEAL::Setting up system... done (0.0213910s)
1619DEAL:: Number of degrees of freedom: 51842
1620DEAL::Assembling system matrix... done (0.0414300s)
1621DEAL::Solving linear system... done (1.56621s)
1622DEAL::Generating output... done (0.729605s)
1623[100%] Built target run
1624
1625> make run
1626[ 66%] Built target @ref step_29 "step-29"
1627[100%] Run @ref step_29 "step-29" with Release configuration
1628DEAL::Generating grid... done (0.0668490s)
1629DEAL:: Number of active cells: 102400
1630DEAL::Setting up system... done (0.109694s)
1631DEAL:: Number of degrees of freedom: 206082
1632DEAL::Assembling system matrix... done (0.160784s)
1633DEAL::Solving linear system... done (7.86577s)
1634DEAL::Generating output... done (2.89320s)
1635[100%] Built target run
1636
1637> make run
1638[ 66%] Built target @ref step_29 "step-29"
1639[100%] Run @ref step_29 "step-29" with Release configuration
1640DEAL::Generating grid... done (0.293154s)
1641DEAL:: Number of active cells: 409600
1642DEAL::Setting up system... done (0.491301s)
1643DEAL:: Number of degrees of freedom: 821762
1644DEAL::Assembling system matrix... done (0.605386s)
1645DEAL::Solving linear system... done (45.1989s)
1646DEAL::Generating output... done (11.2292s)
1647[100%] Built target run
1648@endcode
1649
1650Each time we refine the mesh once, so the number of cells and degrees
1651of freedom roughly quadruples from each step to the next. As can be seen,
1652generating the grid, setting up degrees of freedom, assembling the
1653linear system, and generating output scale pretty closely to linear,
1654whereas solving the linear system is an operation that requires 8
1655times more time each time the number of degrees of freedom is
1656increased by a factor of 4, i.e. it is @f${\cal O}(N^{3/2})@f$. This can
1657be explained by the fact that (using optimal ordering) the
1658bandwidth of a finite element matrix is @f$B={\cal O}(N^{(dim-1)/dim})@f$,
1659and the effort to solve a banded linear system using LU decomposition
1660is @f${\cal O}(BN)@f$. This also explains why the program does run in 3d
1661as well (after changing the dimension on the
1662<code>UltrasoundProblem</code> object), but scales very badly and
1663takes extraordinary patience before it finishes solving the linear
1664system on a mesh with appreciable resolution, even though all the
1665other parts of the program scale very nicely.
1666
1667
1668
1669<a name="extensions"></a>
1670<a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1671
1672
1673An obvious possible extension for this program is to run it in 3d
1674&mdash; after all, the world around us is three-dimensional, and
1675ultrasound beams propagate in three-dimensional media. You can try
1676this by simply changing the template parameter of the principal class
1677in <code>main()</code> and running it. This won't get you very far,
1678though: certainly not if you do 5 global refinement steps as set in
1679the parameter file. You'll simply run out of memory as both the mesh
1680(with its @f$(2^5)^3 \cdot 5^3=2^{15}\cdot 125 \approx 4\cdot 10^6@f$ cells)
1681and in particular the sparse direct solver take too much memory. You
1682can solve with 3 global refinement steps, however, if you have a bit
1683of time: in early 2011, the direct solve takes about half an
1684hour. What you'll notice, however, is that the solution is completely
1685wrong: the mesh size is simply not small enough to resolve the
1686solution's waves accurately, and you can see this in plots of the
1687solution. Consequently, this is one of the cases where adaptivity is
1688indispensable if you don't just want to throw a bigger (presumably
1689%parallel) machine at the problem.
1690 *
1691 *
1692<a name="PlainProg"></a>
1693<h1> The plain program</h1>
1694@include "step-29.cc"
1695*/
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=std::vector< DataComponentInterpretation::DataComponentInterpretation >())
virtual void evaluate_vector_field(const DataPostprocessorInputs::Vector< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
virtual void evaluate_scalar_field(const DataPostprocessorInputs::Scalar< dim > &input_data, std::vector< Vector< double > > &computed_quantities) const
Definition: fe_q.h:549
void enter_subsection(const std::string &subsection)
Definition: point.h:111
Definition: timer.h:119
double cpu_time() const
Definition: timer.cc:236
double stop()
Definition: timer.cc:194
Definition: vector.h:110
VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &x)
Point< 3 > center
UpdateFlags
@ update_hessians
Second derivatives of shape functions.
@ update_values
Shape function values.
@ update_normal_vectors
Normal vectors.
@ update_gradients
Shape function gradients.
Point< 2 > second
Definition: grid_out.cc:4588
Point< 2 > first
Definition: grid_out.cc:4587
unsigned int level
Definition: grid_out.cc:4590
__global__ void set(Number *val, const Number s, const size_type N)
static void declare_parameters(ParameterHandler &prm)
#define Assert(cond, exc)
Definition: exceptions.h:1465
static ::ExceptionBase & ExcDimensionMismatch(std::size_t arg1, std::size_t arg2)
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())
Definition: loop.h:439
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)
const Event initial
Definition: event.cc:65
void subdivided_hyper_cube(Triangulation< dim, spacedim > &tria, const unsigned int repetitions, const double left=0., const double right=1., const bool colorize=false)
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:2042
static const types::blas_int zero
@ matrix
Contents is actually a matrix.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
static const types::blas_int one
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim > > > &Du)
Definition: divergence.h:472
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > w(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
VectorType::value_type * end(VectorType &V)
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)
Definition: work_stream.h:472
int(&) functions(const void *v1, const void *v2)
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition: types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation