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