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