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