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-72.h
Go to the documentation of this file.
1
488) const
489 *   {
490 *   return std::sin(2 * numbers::PI * (p[0] + p[1]));
491 *   }
492 *  
493 *  
494 * @endcode
495 *
496 *
497 * <a name="step_72-ThecodeMinimalSurfaceProblemcodeclassimplementation"></a>
498 * <h3>The <code>MinimalSurfaceProblem</code> class implementation</h3>
499 *
500
501 *
502 *
503 * <a name="step_72-MinimalSurfaceProblemMinimalSurfaceProblem"></a>
504 * <h4>MinimalSurfaceProblem::MinimalSurfaceProblem</h4>
505 *
506
507 *
508 * There have been no changes made to the class constructor.
509 *
510 * @code
511 *   template <int dim>
512 *   MinimalSurfaceProblem<dim>::MinimalSurfaceProblem()
513 *   : dof_handler(triangulation)
514 *   , fe(2)
515 *   , quadrature_formula(fe.degree + 1)
516 *   {}
517 *  
518 *  
519 * @endcode
520 *
521 *
522 * <a name="step_72-MinimalSurfaceProblemsetup_system"></a>
523 * <h4>MinimalSurfaceProblem::setup_system</h4>
524 *
525
526 *
527 * There have been no changes made to the function that sets up the class
528 * data structures, namely the DoFHandler, the hanging node constraints
529 * applied to the problem, and the linear system.
530 *
531 * @code
532 *   template <int dim>
533 *   void MinimalSurfaceProblem<dim>::setup_system()
534 *   {
535 *   dof_handler.distribute_dofs(fe);
536 *   current_solution.reinit(dof_handler.n_dofs());
537 *   newton_update.reinit(dof_handler.n_dofs());
538 *   system_rhs.reinit(dof_handler.n_dofs());
539 *  
540 *   zero_constraints.clear();
542 *   0,
544 *   zero_constraints);
545 *   DoFTools::make_hanging_node_constraints(dof_handler, zero_constraints);
546 *   zero_constraints.close();
547 *  
548 *   nonzero_constraints.clear();
550 *   0,
551 *   BoundaryValues<dim>(),
552 *   nonzero_constraints);
553 *   nonzero_constraints.close();
554 *  
555 *   DynamicSparsityPattern dsp(dof_handler.n_dofs());
556 *   DoFTools::make_sparsity_pattern(dof_handler, dsp, zero_constraints);
557 *  
558 *   sparsity_pattern.copy_from(dsp);
559 *   system_matrix.reinit(sparsity_pattern);
560 *   }
561 *  
562 * @endcode
563 *
564 *
565 * <a name="step_72-Assemblingthelinearsystem"></a>
566 * <h4>Assembling the linear system</h4>
567 *
568
569 *
570 *
571 * <a name="step_72-Manualassembly"></a>
572 * <h5>Manual assembly</h5>
573 *
574
575 *
576 * The assembly functions are the interesting contributions to this tutorial.
577 * The assemble_system_unassisted() method implements exactly the same
578 * assembly function as is detailed in @ref step_15 "step-15", but in this instance we
579 * use the MeshWorker::mesh_loop() function to multithread the assembly
580 * process. The reason for doing this is quite simple: When using
581 * automatic differentiation, we know that there is to be some additional
582 * computational overhead incurred. In order to mitigate this performance
583 * loss, we'd like to take advantage of as many (easily available)
584 * computational resources as possible. The MeshWorker::mesh_loop() concept
585 * makes this a relatively straightforward task. At the same time, for the
586 * purposes of fair comparison, we need to do the same to the implementation
587 * that uses no assistance when computing the residual or its linearization.
588 * (The MeshWorker::mesh_loop() function is first discussed in @ref step_12 "step-12" and
589 * @ref step_16 "step-16", if you'd like to read up on it.)
590 *
591
592 *
593 * The steps required to implement the multithreading are the same between the
594 * three functions, so we'll use the assemble_system_unassisted() function
595 * as an opportunity to focus on the multithreading itself.
596 *
597 * @code
598 *   template <int dim>
599 *   void MinimalSurfaceProblem<dim>::assemble_system_unassisted()
600 *   {
601 *   system_matrix = 0;
602 *   system_rhs = 0;
603 *  
604 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
605 *  
606 * @endcode
607 *
608 * The MeshWorker::mesh_loop() expects that we provide two exemplar data
609 * structures. The first, `ScratchData`, is to store all large data that
610 * is to be reused between threads. The `CopyData` will hold the
611 * contributions to the linear system that come from each cell. These
612 * independent matrix-vector pairs must be accumulated into the
613 * global linear system sequentially. Since we don't need anything
614 * on top of what the MeshWorker::ScratchData and MeshWorker::CopyData
615 * classes already provide, we use these exact class definitions for
616 * our problem. Note that we only require a single instance of a local
617 * matrix, local right-hand side vector, and cell degree of freedom index
618 * vector -- the MeshWorker::CopyData therefore has `1` for all three
619 * of its template arguments.
620 *
621 * @code
622 *   using ScratchData = MeshWorker::ScratchData<dim>;
623 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
624 *  
625 * @endcode
626 *
627 * We also need to know what type of iterator we'll be working with
628 * during assembly. For simplicity, we just ask the compiler to work
629 * this out for us using the decltype() specifier, knowing that we'll
630 * be iterating over active cells owned by the @p dof_handler.
631 *
632 * @code
633 *   using CellIteratorType = decltype(dof_handler.begin_active());
634 *  
635 * @endcode
636 *
637 * Here we initialize the exemplar data structures. Since we know that
638 * we need to compute the shape function gradients, weighted Jacobian,
639 * and the position of the quadrate points in real space, we pass these
640 * flags into the class constructor.
641 *
642 * @code
643 *   const ScratchData sample_scratch_data(fe,
644 *   quadrature_formula,
645 *   update_gradients |
647 *   update_JxW_values);
648 *   const CopyData sample_copy_data(dofs_per_cell);
649 *  
650 * @endcode
651 *
652 * Now we define a lambda function that will perform the assembly on
653 * a single cell. The three arguments are those that will be expected by
654 * MeshWorker::mesh_loop(), due to the arguments that we'll pass to that
655 * final call. We also capture the @p this pointer, which means that we'll
656 * have access to "this" (i.e., the current `MinimalSurfaceProblem<dim>`)
657 * class instance, and its private member data (since the lambda function is
658 * defined within a MinimalSurfaceProblem<dim> method).
659 *
660
661 *
662 * At the top of the function, we initialize the data structures
663 * that are dependent on the cell for which the work is being
664 * performed. Observe that the reinitialization call actually
665 * returns an instance to an FEValues object that is initialized
666 * and stored within (and, therefore, reused by) the
667 * `scratch_data` object.
668 *
669
670 *
671 * Similarly, we get aliases to the local matrix, local RHS
672 * vector, and local cell DoF indices from the `copy_data`
673 * instance that MeshWorker::mesh_loop() provides. We then
674 * initialize the cell DoF indices, knowing that the local matrix
675 * and vector are already correctly sized.
676 *
677 * @code
678 *   const auto cell_worker = [this](const CellIteratorType &cell,
679 *   ScratchData &scratch_data,
680 *   CopyData &copy_data) {
681 *   const auto &fe_values = scratch_data.reinit(cell);
682 *  
683 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
684 *   Vector<double> &cell_rhs = copy_data.vectors[0];
685 *   std::vector<types::global_dof_index> &local_dof_indices =
686 *   copy_data.local_dof_indices[0];
687 *   cell->get_dof_indices(local_dof_indices);
688 *  
689 * @endcode
690 *
691 * For Newton's method, we require the gradient of the solution at the
692 * point about which the problem is being linearized.
693 *
694
695 *
696 * Once we have that, we can perform assembly for this cell in
697 * the usual way. One minor difference to @ref step_15 "step-15" is that we've
698 * used the (rather convenient) range-based loops to iterate
699 * over all quadrature points and degrees-of-freedom.
700 *
701 * @code
702 *   std::vector<Tensor<1, dim>> old_solution_gradients(
703 *   fe_values.n_quadrature_points);
704 *   fe_values.get_function_gradients(current_solution,
705 *   old_solution_gradients);
706 *  
707 *   for (const unsigned int q : fe_values.quadrature_point_indices())
708 *   {
709 *   const double coeff =
710 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
711 *   old_solution_gradients[q]);
712 *  
713 *   for (const unsigned int i : fe_values.dof_indices())
714 *   {
715 *   for (const unsigned int j : fe_values.dof_indices())
716 *   cell_matrix(i, j) +=
717 *   (((fe_values.shape_grad(i, q) // ((\nabla \phi_i
718 *   * coeff // * a_n
719 *   * fe_values.shape_grad(j, q)) // * \nabla \phi_j)
720 *   - // -
721 *   (fe_values.shape_grad(i, q) // (\nabla \phi_i
722 *   * coeff * coeff * coeff // * a_n^3
723 *   * (fe_values.shape_grad(j, q) // * (\nabla \phi_j
724 *   * old_solution_gradients[q]) // * \nabla u_n)
725 *   * old_solution_gradients[q])) // * \nabla u_n)))
726 *   * fe_values.JxW(q)); // * dx
727 *  
728 *   cell_rhs(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
729 *   * coeff // * a_n
730 *   * old_solution_gradients[q] // * \nabla u_n
731 *   * fe_values.JxW(q)); // * dx
732 *   }
733 *   }
734 *   };
735 *  
736 * @endcode
737 *
738 * The second lambda function that MeshWorker::mesh_loop() requires is
739 * one that performs the task of accumulating the local contributions
740 * in the global linear system. That is precisely what this one does,
741 * and the details of the implementation have been seen before. The
742 * primary point to recognize is that the local contributions are stored
743 * in the `copy_data` instance that is passed into this function. This
744 * `copy_data` has been filled with data during @a some call to the
745 * `cell_worker`.
746 *
747 * @code
748 *   const auto copier = [this](const CopyData &copy_data) {
749 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
750 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
751 *   const std::vector<types::global_dof_index> &local_dof_indices =
752 *   copy_data.local_dof_indices[0];
753 *  
754 *   zero_constraints.distribute_local_to_global(
755 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
756 *   };
757 *  
758 * @endcode
759 *
760 * We have all of the required functions definitions in place, so
761 * now we call the MeshWorker::mesh_loop() to perform the actual
762 * assembly. We pass a flag as the last parameter which states
763 * that we only want to perform the assembly on the
764 * cells. Internally, MeshWorker::mesh_loop() then distributes the
765 * available work to different threads, making efficient use of
766 * the multiple cores almost all of today's processors have to
767 * offer.
768 *
769 * @code
770 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
771 *   cell_worker,
772 *   copier,
773 *   sample_scratch_data,
774 *   sample_copy_data,
775 *   MeshWorker::assemble_own_cells);
776 *   }
777 *  
778 * @endcode
779 *
780 *
781 * <a name="step_72-Assemblyviadifferentiationoftheresidualvector"></a>
782 * <h5>Assembly via differentiation of the residual vector</h5>
783 *
784
785 *
786 * As outlined in the introduction, what we need to do for this
787 * second approach is implement the local contributions @f$F(U)^K@f$
788 * from cell @f$K@f$ to the residual vector, and then let the
789 * AD machinery deal with how to compute the
790 * derivatives @f$J(U)_{ij}^K=\frac{\partial F(U)^K_i}{\partial U_j}@f$
791 * from it.
792 *
793
794 *
795 * For the following, recall that
796 * @f[
797 * F(U)_i^K \dealcoloneq
798 * \int\limits_K\nabla \varphi_i \cdot \left[ \frac{1}{\sqrt{1+|\nabla
799 * u|^{2}}} \nabla u \right] \, dV ,
800 * @f]
801 * where @f$u(\mathbf x)=\sum_j U_j \varphi_j(\mathbf x)@f$.
802 *
803
804 *
805 * Let us see how this is implemented in practice:
806 *
807 * @code
808 *   template <int dim>
809 *   void MinimalSurfaceProblem<dim>::assemble_system_with_residual_linearization()
810 *   {
811 *   system_matrix = 0;
812 *   system_rhs = 0;
813 *  
814 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
815 *  
816 *   using ScratchData = MeshWorker::ScratchData<dim>;
817 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
818 *   using CellIteratorType = decltype(dof_handler.begin_active());
819 *  
820 *   const ScratchData sample_scratch_data(fe,
821 *   quadrature_formula,
822 *   update_gradients |
824 *   update_JxW_values);
825 *   const CopyData sample_copy_data(dofs_per_cell);
826 *  
827 * @endcode
828 *
829 * We'll define up front the AD data structures that we'll be using,
830 * utilizing the techniques shown in @ref step_71 "step-71".
831 * In this case, we choose the helper class that will automatically compute
832 * the linearization of the finite element residual using Sacado forward
833 * automatic differentiation types. These number types can be used to
834 * compute first derivatives only. This is exactly what we want, because we
835 * know that we'll only be linearizing the residual, which means that we
836 * only need to compute first-order derivatives. The return values from the
837 * calculations are to be of type `double`.
838 *
839
840 *
841 * We also need an extractor to retrieve some data related to the field
842 * solution to the problem.
843 *
844 * @code
845 *   using ADHelper = Differentiation::AD::ResidualLinearization<
846 *   Differentiation::AD::NumberTypes::sacado_dfad,
847 *   double>;
848 *   using ADNumberType = typename ADHelper::ad_type;
849 *  
850 *   const FEValuesExtractors::Scalar u_fe(0);
851 *  
852 * @endcode
853 *
854 * With this, let us define the lambda function that will be used
855 * to compute the cell contributions to the Jacobian matrix and
856 * the right hand side:
857 *
858 * @code
859 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
860 *   ScratchData &scratch_data,
861 *   CopyData &copy_data) {
862 *   const auto &fe_values = scratch_data.reinit(cell);
863 *   const unsigned int dofs_per_cell = fe_values.get_fe().n_dofs_per_cell();
864 *  
865 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
866 *   Vector<double> &cell_rhs = copy_data.vectors[0];
867 *   std::vector<types::global_dof_index> &local_dof_indices =
868 *   copy_data.local_dof_indices[0];
869 *   cell->get_dof_indices(local_dof_indices);
870 *  
871 * @endcode
872 *
873 * We'll now create and initialize an instance of the AD helper class.
874 * To do this, we need to specify how many independent variables and
875 * dependent variables there are. The independent variables will be the
876 * number of local degrees of freedom that our solution vector has,
877 * i.e., the number @f$j@f$ in the per-element representation of the
878 * discretized solution vector
879 * @f$u (\mathbf{x})|_K = \sum\limits_{j} U^K_i \varphi_j(\mathbf{x})@f$
880 * that indicates how many solution coefficients are associated with
881 * each finite element. In deal.II, this equals
882 * FiniteElement::dofs_per_cell. The number of dependent variables will be
883 * the number of entries in the local residual vector that we will be
884 * forming. In this particular problem (like many others that employ the
885 * [standard Galerkin
886 * method](https://en.wikipedia.org/wiki/Galerkin_method)) the number of
887 * local solution coefficients matches the number of local residual
888 * equations.
889 *
890 * @code
891 *   const unsigned int n_independent_variables = local_dof_indices.size();
892 *   const unsigned int n_dependent_variables = dofs_per_cell;
893 *   ADHelper ad_helper(n_independent_variables, n_dependent_variables);
894 *  
895 * @endcode
896 *
897 * Next we inform the helper of the values of the solution, i.e., the
898 * actual values for @f$U_j@f$ about which we
899 * wish to linearize. As this is done on each element individually, we
900 * have to extract the solution coefficients from the global solution
901 * vector. In other words, we define all of those coefficients @f$U_j@f$
902 * where @f$j@f$ is a local degree of freedom as the independent variables
903 * that enter the computation of the vector @f$F(U)^{K}@f$ (the dependent
904 * function).
905 *
906
907 *
908 * Then we get the complete set of degree of freedom values as
909 * represented by auto-differentiable numbers. The operations
910 * performed with these variables are tracked by the AD library
911 * from this point until the object goes out of scope. So it is
912 * <em>precisely these variables</em> with respect to which we will
913 * compute derivatives of the residual entries.
914 *
915 * @code
916 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
917 *  
918 *   const std::vector<ADNumberType> &dof_values_ad =
919 *   ad_helper.get_sensitive_dof_values();
920 *  
921 * @endcode
922 *
923 * Then we do some problem specific tasks, the first being to
924 * compute all values, (spatial) gradients, and the like based on
925 * "sensitive" AD degree of freedom values. In this instance we are
926 * retrieving the solution gradients at each quadrature point. Observe
927 * that the solution gradients are now sensitive
928 * to the values of the degrees of freedom as they use the @p ADNumberType
929 * as the scalar type and the @p dof_values_ad vector provides the local
930 * DoF values.
931 *
932 * @code
933 *   std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
934 *   fe_values.n_quadrature_points);
935 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
936 *   dof_values_ad, old_solution_gradients);
937 *  
938 * @endcode
939 *
940 * The next variable that we declare will store the cell residual vector
941 * contributions. This is rather self-explanatory, save for one
942 * <b>very important</b> detail:
943 * Note that each entry in the vector is hand-initialized with a value
944 * of zero. This is a <em>highly recommended</em> practice, as some AD
945 * libraries appear not to safely initialize the internal data
946 * structures of these number types. Not doing so could lead to some
947 * very hard to understand or detect bugs (appreciate that the author
948 * of this program mentions this out of, generally bad, experience). So
949 * out of an abundance of caution it's worthwhile zeroing the initial
950 * value explicitly. After that, apart from a sign change the residual
951 * assembly looks much the same as we saw for the cell RHS vector before:
952 * We loop over all quadrature points, ensure that the coefficient now
953 * encodes its dependence on the (sensitive) finite element DoF values by
954 * using the correct `ADNumberType`, and finally we assemble the
955 * components of the residual vector. For complete clarity, the finite
956 * element shape functions (and their gradients, etc.) as well as the
957 * "JxW" values remain scalar
958 * valued, but the @p coeff and the @p old_solution_gradients at each
959 * quadrature point are computed in terms of the independent
960 * variables.
961 *
962 * @code
963 *   std::vector<ADNumberType> residual_ad(n_dependent_variables,
964 *   ADNumberType(0.0));
965 *   for (const unsigned int q : fe_values.quadrature_point_indices())
966 *   {
967 *   const ADNumberType coeff =
968 *   1.0 / std::sqrt(1.0 + old_solution_gradients[q] *
969 *   old_solution_gradients[q]);
970 *  
971 *   for (const unsigned int i : fe_values.dof_indices())
972 *   {
973 *   residual_ad[i] += (fe_values.shape_grad(i, q) // \nabla \phi_i
974 *   * coeff // * a_n
975 *   * old_solution_gradients[q]) // * \nabla u_n
976 *   * fe_values.JxW(q); // * dx
977 *   }
978 *   }
979 *  
980 * @endcode
981 *
982 * Once we have the full cell residual vector computed, we can register
983 * it with the helper class.
984 *
985
986 *
987 * Thereafter, we compute the residual values (basically,
988 * extracting the real values from what we already computed) and
989 * their Jacobian (the linearization of each residual component
990 * with respect to all cell DoFs) at the evaluation point. For
991 * the purposes of assembly into the global linear system, we
992 * have to respect the sign difference between the residual and
993 * the RHS contribution: For Newton's method, the right hand
994 * side vector needs to be equal to the *negative* residual
995 * vector.
996 *
997 * @code
998 *   ad_helper.register_residual_vector(residual_ad);
999 *  
1000 *   ad_helper.compute_residual(cell_rhs);
1001 *   cell_rhs *= -1.0;
1002 *  
1003 *   ad_helper.compute_linearization(cell_matrix);
1004 *   };
1005 *  
1006 * @endcode
1007 *
1008 * The remainder of the function equals what we had previously:
1009 *
1010 * @code
1011 *   const auto copier = [this](const CopyData &copy_data) {
1012 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1013 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
1014 *   const std::vector<types::global_dof_index> &local_dof_indices =
1015 *   copy_data.local_dof_indices[0];
1016 *  
1017 *   zero_constraints.distribute_local_to_global(
1018 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1019 *   };
1020 *  
1021 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1022 *   cell_worker,
1023 *   copier,
1024 *   sample_scratch_data,
1025 *   sample_copy_data,
1027 *   }
1028 *  
1029 * @endcode
1030 *
1031 *
1032 * <a name="step_72-Assemblyviadifferentiationoftheenergyfunctional"></a>
1033 * <h5>Assembly via differentiation of the energy functional</h5>
1034 *
1035
1036 *
1037 * In this third approach, we compute residual and Jacobian as first
1038 * and second derivatives of the local energy functional
1039 * @f[
1040 * E\left( U \right)^K
1041 * \dealcoloneq \int\limits_{K} \Psi \left( u \right) \, dV
1042 * \approx \sum\limits_{q}^{n_{\textrm{q-points}}} \Psi \left( u \left(
1043 * \mathbf{X}_{q} \right) \right) \underbrace{\vert J_{q} \vert \times
1044 * W_{q}}_{\text{JxW(q)}}
1045 * @f]
1046 * with the energy density given by
1047 * @f[
1048 * \Psi \left( u \right) = \sqrt{1+|\nabla u|^{2}} .
1049 * @f]
1050 *
1051
1052 *
1053 * Let us again see how this is done:
1054 *
1055 * @code
1056 *   template <int dim>
1057 *   void MinimalSurfaceProblem<dim>::assemble_system_using_energy_functional()
1058 *   {
1059 *   system_matrix = 0;
1060 *   system_rhs = 0;
1061 *  
1062 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1063 *  
1064 *   using ScratchData = MeshWorker::ScratchData<dim>;
1065 *   using CopyData = MeshWorker::CopyData<1, 1, 1>;
1066 *   using CellIteratorType = decltype(dof_handler.begin_active());
1067 *  
1068 *   const ScratchData sample_scratch_data(fe,
1069 *   quadrature_formula,
1070 *   update_gradients |
1072 *   update_JxW_values);
1073 *   const CopyData sample_copy_data(dofs_per_cell);
1074 *  
1075 * @endcode
1076 *
1077 * In this implementation of the assembly process, we choose the helper
1078 * class that will automatically compute both the residual and its
1079 * linearization from the cell contribution to an energy functional using
1080 * nested Sacado forward automatic differentiation types.
1081 * The selected number types can be used to compute both first and
1082 * second derivatives. We require this, as the residual defined as the
1083 * sensitivity of the potential energy with respect to the DoF values (i.e.
1084 * its gradient). We'll then need to linearize the residual, implying that
1085 * second derivatives of the potential energy must be computed. You might
1086 * want to compare this with the definition of `ADHelper` used int
1087 * previous function, where we used
1088 * `Differentiation::AD::ResidualLinearization<Differentiation::AD::NumberTypes::sacado_dfad,double>`.
1089 *
1090 * @code
1091 *   using ADHelper = Differentiation::AD::EnergyFunctional<
1092 *   Differentiation::AD::NumberTypes::sacado_dfad_dfad,
1093 *   double>;
1094 *   using ADNumberType = typename ADHelper::ad_type;
1095 *  
1096 *   const FEValuesExtractors::Scalar u_fe(0);
1097 *  
1098 * @endcode
1099 *
1100 * Let us then again define the lambda function that does the integration on
1101 * a cell.
1102 *
1103
1104 *
1105 * To initialize an instance of the helper class, we now only require
1106 * that the number of independent variables (that is, the number
1107 * of degrees of freedom associated with the element solution vector)
1108 * are known up front. This is because the second-derivative matrix that
1109 * results from an energy functional is necessarily square (and also,
1110 * incidentally, symmetric).
1111 *
1112 * @code
1113 *   const auto cell_worker = [&u_fe, this](const CellIteratorType &cell,
1114 *   ScratchData &scratch_data,
1115 *   CopyData &copy_data) {
1116 *   const auto &fe_values = scratch_data.reinit(cell);
1117 *  
1118 *   FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1119 *   Vector<double> &cell_rhs = copy_data.vectors[0];
1120 *   std::vector<types::global_dof_index> &local_dof_indices =
1121 *   copy_data.local_dof_indices[0];
1122 *   cell->get_dof_indices(local_dof_indices);
1123 *  
1124 *   const unsigned int n_independent_variables = local_dof_indices.size();
1125 *   ADHelper ad_helper(n_independent_variables);
1126 *  
1127 * @endcode
1128 *
1129 * Once more, we register all cell DoFs values with the helper, followed
1130 * by extracting the "sensitive" variant of these values that are to be
1131 * used in subsequent operations that must be differentiated -- one of
1132 * those being the calculation of the solution gradients.
1133 *
1134 * @code
1135 *   ad_helper.register_dof_values(current_solution, local_dof_indices);
1136 *  
1137 *   const std::vector<ADNumberType> &dof_values_ad =
1138 *   ad_helper.get_sensitive_dof_values();
1139 *  
1140 *   std::vector<Tensor<1, dim, ADNumberType>> old_solution_gradients(
1141 *   fe_values.n_quadrature_points);
1142 *   fe_values[u_fe].get_function_gradients_from_local_dof_values(
1143 *   dof_values_ad, old_solution_gradients);
1144 *  
1145 * @endcode
1146 *
1147 * We next create a variable that stores the cell total energy.
1148 * Once more we emphasize that we explicitly zero-initialize this value,
1149 * thereby ensuring the integrity of the data for this starting value.
1150 *
1151
1152 *
1153 * The aim for our approach is then to compute the cell total
1154 * energy, which is the sum of the internal (due to right hand
1155 * side functions, typically linear in @f$U@f$) and external
1156 * energies. In this particular case, we have no external
1157 * energies (e.g., from source terms or Neumann boundary
1158 * conditions), so we'll focus on the internal energy part.
1159 *
1160
1161 *
1162 * In fact, computing @f$E(U)^K@f$ is almost trivial, requiring only
1163 * the following lines:
1164 *
1165 * @code
1166 *   ADNumberType energy_ad = ADNumberType(0.0);
1167 *   for (const unsigned int q : fe_values.quadrature_point_indices())
1168 *   {
1169 *   const ADNumberType psi = std::sqrt(1.0 + old_solution_gradients[q] *
1170 *   old_solution_gradients[q]);
1171 *  
1172 *   energy_ad += psi * fe_values.JxW(q);
1173 *   }
1174 *  
1175 * @endcode
1176 *
1177 * After we've computed the total energy on this cell, we'll
1178 * register it with the helper. Based on that, we may now
1179 * compute the desired quantities, namely the residual values
1180 * and their Jacobian at the evaluation point. As before, the
1181 * Newton right hand side needs to be the negative of the
1182 * residual:
1183 *
1184 * @code
1185 *   ad_helper.register_energy_functional(energy_ad);
1186 *  
1187 *   ad_helper.compute_residual(cell_rhs);
1188 *   cell_rhs *= -1.0;
1189 *  
1190 *   ad_helper.compute_linearization(cell_matrix);
1191 *   };
1192 *  
1193 * @endcode
1194 *
1195 * As in the previous two functions, the remainder of the function is as
1196 * before:
1197 *
1198 * @code
1199 *   const auto copier = [this](const CopyData &copy_data) {
1200 *   const FullMatrix<double> &cell_matrix = copy_data.matrices[0];
1201 *   const Vector<double> &cell_rhs = copy_data.vectors[0];
1202 *   const std::vector<types::global_dof_index> &local_dof_indices =
1203 *   copy_data.local_dof_indices[0];
1204 *  
1205 *   zero_constraints.distribute_local_to_global(
1206 *   cell_matrix, cell_rhs, local_dof_indices, system_matrix, system_rhs);
1207 *   };
1208 *  
1209 *   MeshWorker::mesh_loop(dof_handler.active_cell_iterators(),
1210 *   cell_worker,
1211 *   copier,
1212 *   sample_scratch_data,
1213 *   sample_copy_data,
1215 *   }
1216 *  
1217 *  
1218 * @endcode
1219 *
1220 *
1221 * <a name="step_72-MinimalSurfaceProblemsolve"></a>
1222 * <h4>MinimalSurfaceProblem::solve</h4>
1223 *
1224
1225 *
1226 * The solve function is the same as is used in @ref step_15 "step-15".
1227 *
1228 * @code
1229 *   template <int dim>
1230 *   void MinimalSurfaceProblem<dim>::solve()
1231 *   {
1232 *   SolverControl solver_control(system_rhs.size(),
1233 *   system_rhs.l2_norm() * 1e-6);
1234 *   SolverCG<Vector<double>> solver(solver_control);
1235 *  
1236 *   PreconditionSSOR<SparseMatrix<double>> preconditioner;
1237 *   preconditioner.initialize(system_matrix, 1.2);
1238 *  
1239 *   solver.solve(system_matrix, newton_update, system_rhs, preconditioner);
1240 *  
1241 *   zero_constraints.distribute(newton_update);
1242 *  
1243 *   const double alpha = determine_step_length();
1244 *   current_solution.add(alpha, newton_update);
1245 *   }
1246 *  
1247 *  
1248 * @endcode
1249 *
1250 *
1251 * <a name="step_72-MinimalSurfaceProblemrefine_mesh"></a>
1252 * <h4>MinimalSurfaceProblem::refine_mesh</h4>
1253 *
1254
1255 *
1256 * Nothing has changed since @ref step_15 "step-15" with respect to the mesh refinement
1257 * procedure and transfer of the solution between adapted meshes.
1258 *
1259 * @code
1260 *   template <int dim>
1261 *   void MinimalSurfaceProblem<dim>::refine_mesh()
1262 *   {
1263 *   Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1264 *  
1266 *   dof_handler,
1267 *   QGauss<dim - 1>(fe.degree + 1),
1268 *   std::map<types::boundary_id, const Function<dim> *>(),
1269 *   current_solution,
1270 *   estimated_error_per_cell);
1271 *  
1273 *   estimated_error_per_cell,
1274 *   0.3,
1275 *   0.03);
1276 *  
1278 *  
1279 *   SolutionTransfer<dim> solution_transfer(dof_handler);
1280 *   const Vector<double> coarse_solution = current_solution;
1281 *   solution_transfer.prepare_for_coarsening_and_refinement(coarse_solution);
1282 *  
1284 *  
1285 *   setup_system();
1286 *  
1287 *   solution_transfer.interpolate(coarse_solution, current_solution);
1288 *   nonzero_constraints.distribute(current_solution);
1289 *   }
1290 *  
1291 *  
1292 *  
1293 * @endcode
1294 *
1295 *
1296 * <a name="step_72-MinimalSurfaceProblemcompute_residual"></a>
1297 * <h4>MinimalSurfaceProblem::compute_residual</h4>
1298 *
1299
1300 *
1301 * ... as does the function used to compute the residual during the
1302 * solution iteration procedure. One could replace this by
1303 * differentiation of the energy functional if one really wanted,
1304 * but for simplicity we here simply copy what we already had in
1305 * @ref step_15 "step-15".
1306 *
1307 * @code
1308 *   template <int dim>
1309 *   double MinimalSurfaceProblem<dim>::compute_residual(const double alpha) const
1310 *   {
1311 *   Vector<double> residual(dof_handler.n_dofs());
1312 *  
1313 *   Vector<double> evaluation_point(dof_handler.n_dofs());
1314 *   evaluation_point = current_solution;
1315 *   evaluation_point.add(alpha, newton_update);
1316 *  
1317 *   FEValues<dim> fe_values(fe,
1318 *   quadrature_formula,
1320 *   update_JxW_values);
1321 *  
1322 *   const unsigned int dofs_per_cell = fe.n_dofs_per_cell();
1323 *   const unsigned int n_q_points = quadrature_formula.size();
1324 *  
1325 *   Vector<double> cell_residual(dofs_per_cell);
1326 *   std::vector<Tensor<1, dim>> gradients(n_q_points);
1327 *  
1328 *   std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
1329 *  
1330 *   for (const auto &cell : dof_handler.active_cell_iterators())
1331 *   {
1332 *   cell_residual = 0;
1333 *   fe_values.reinit(cell);
1334 *  
1335 *   fe_values.get_function_gradients(evaluation_point, gradients);
1336 *  
1337 *   for (unsigned int q = 0; q < n_q_points; ++q)
1338 *   {
1339 *   const double coeff =
1340 *   1.0 / std::sqrt(1.0 + gradients[q] * gradients[q]);
1341 *  
1342 *   for (unsigned int i = 0; i < dofs_per_cell; ++i)
1343 *   cell_residual(i) -= (fe_values.shape_grad(i, q) // \nabla \phi_i
1344 *   * coeff // * a_n
1345 *   * gradients[q] // * \nabla u_n
1346 *   * fe_values.JxW(q)); // * dx
1347 *   }
1348 *  
1349 *   cell->get_dof_indices(local_dof_indices);
1350 *   zero_constraints.distribute_local_to_global(cell_residual,
1351 *   local_dof_indices,
1352 *   residual);
1353 *   }
1354 *  
1355 *   return residual.l2_norm();
1356 *   }
1357 *  
1358 *  
1359 *  
1360 * @endcode
1361 *
1362 *
1363 * <a name="step_72-MinimalSurfaceProblemdetermine_step_length"></a>
1364 * <h4>MinimalSurfaceProblem::determine_step_length</h4>
1365 *
1366
1367 *
1368 * The choice of step length (or, under-relaxation factor) for the nonlinear
1369 * iterations procedure remains fixed at the value chosen and discussed in
1370 * @ref step_15 "step-15".
1371 *
1372 * @code
1373 *   template <int dim>
1374 *   double MinimalSurfaceProblem<dim>::determine_step_length() const
1375 *   {
1376 *   return 0.1;
1377 *   }
1378 *  
1379 *  
1380 *  
1381 * @endcode
1382 *
1383 *
1384 * <a name="step_72-MinimalSurfaceProblemoutput_results"></a>
1385 * <h4>MinimalSurfaceProblem::output_results</h4>
1386 *
1387
1388 *
1389 * This last function to be called from `run()` outputs the current solution
1390 * (and the Newton update) in graphical form as a VTU file. It is entirely the
1391 * same as what has been used in previous tutorials.
1392 *
1393 * @code
1394 *   template <int dim>
1395 *   void MinimalSurfaceProblem<dim>::output_results(
1396 *   const unsigned int refinement_cycle) const
1397 *   {
1398 *   DataOut<dim> data_out;
1399 *  
1400 *   data_out.attach_dof_handler(dof_handler);
1401 *   data_out.add_data_vector(current_solution, "solution");
1402 *   data_out.add_data_vector(newton_update, "update");
1403 *   data_out.build_patches();
1404 *  
1405 *   const std::string filename =
1406 *   "solution-" + Utilities::int_to_string(refinement_cycle, 2) + ".vtu";
1407 *   std::ofstream output(filename);
1408 *   data_out.write_vtu(output);
1409 *   }
1410 *  
1411 *  
1412 * @endcode
1413 *
1414 *
1415 * <a name="step_72-MinimalSurfaceProblemrun"></a>
1416 * <h4>MinimalSurfaceProblem::run</h4>
1417 *
1418
1419 *
1420 * In the run function, most remains the same as was first implemented
1421 * in @ref step_15 "step-15". The only observable changes are that we can now choose (via
1422 * the parameter file) what the final acceptable tolerance for the system
1423 * residual is, and that we can choose which method of assembly we wish to
1424 * utilize. To make the second choice clear, we output to the console some
1425 * message which indicates the selection. Since we're interested in comparing
1426 * the time taken to assemble for each of the three methods, we've also
1427 * added a timer that keeps a track of how much time is spent during assembly.
1428 * We also track the time taken to solve the linear system, so that we can
1429 * contrast those numbers to the part of the code which would normally take
1430 * the longest time to execute.
1431 *
1432 * @code
1433 *   template <int dim>
1434 *   void MinimalSurfaceProblem<dim>::run(const int formulation,
1435 *   const double tolerance)
1436 *   {
1437 *   std::cout << "******** Assembly approach ********" << std::endl;
1438 *   const std::array<std::string, 3> method_descriptions = {
1439 *   {"Unassisted implementation (full hand linearization).",
1440 *   "Automated linearization of the finite element residual.",
1441 *   "Automated computation of finite element residual and linearization using a variational formulation."}};
1442 *   AssertIndexRange(formulation, method_descriptions.size());
1443 *   std::cout << method_descriptions[formulation] << std::endl << std::endl;
1444 *  
1445 *  
1447 *  
1450 *  
1451 *   setup_system();
1452 *   nonzero_constraints.distribute(current_solution);
1453 *  
1454 *   double last_residual_norm = std::numeric_limits<double>::max();
1455 *   unsigned int refinement_cycle = 0;
1456 *   do
1457 *   {
1458 *   std::cout << "Mesh refinement step " << refinement_cycle << std::endl;
1459 *  
1460 *   if (refinement_cycle != 0)
1461 *   refine_mesh();
1462 *  
1463 *   std::cout << " Initial residual: " << compute_residual(0) << std::endl;
1464 *  
1465 *   for (unsigned int inner_iteration = 0; inner_iteration < 5;
1466 *   ++inner_iteration)
1467 *   {
1468 *   {
1469 *   TimerOutput::Scope t(timer, "Assemble");
1470 *  
1471 *   if (formulation == 0)
1472 *   assemble_system_unassisted();
1473 *   else if (formulation == 1)
1474 *   assemble_system_with_residual_linearization();
1475 *   else if (formulation == 2)
1476 *   assemble_system_using_energy_functional();
1477 *   else
1478 *   AssertThrow(false, ExcNotImplemented());
1479 *   }
1480 *  
1481 *   last_residual_norm = system_rhs.l2_norm();
1482 *  
1483 *   {
1484 *   TimerOutput::Scope t(timer, "Solve");
1485 *   solve();
1486 *   }
1487 *  
1488 *  
1489 *   std::cout << " Residual: " << compute_residual(0) << std::endl;
1490 *   }
1491 *  
1492 *   output_results(refinement_cycle);
1493 *  
1494 *   ++refinement_cycle;
1495 *   std::cout << std::endl;
1496 *   }
1497 *   while (last_residual_norm > tolerance);
1498 *   }
1499 *   } // namespace Step72
1500 *  
1501 * @endcode
1502 *
1503 *
1504 * <a name="step_72-Themainfunction"></a>
1505 * <h4>The main function</h4>
1506 *
1507
1508 *
1509 * Finally the main function. This follows the scheme of most other main
1510 * functions, with two obvious exceptions:
1511 * - We call Utilities::MPI::MPI_InitFinalize in order to set up (via a hidden
1512 * default parameter) the number of threads using the execution of
1513 * multithreaded tasks.
1514 * - We also have a few lines dedicates to reading in or initializing the
1515 * user-defined parameters that will be considered during the execution of the
1516 * program.
1517 *
1518 * @code
1519 *   int main(int argc, char *argv[])
1520 *   {
1521 *   try
1522 *   {
1523 *   using namespace Step72;
1524 *  
1525 *   Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);
1526 *  
1527 *   std::string prm_file;
1528 *   if (argc > 1)
1529 *   prm_file = argv[1];
1530 *   else
1531 *   prm_file = "parameters.prm";
1532 *  
1533 *   const MinimalSurfaceProblemParameters parameters;
1534 *   ParameterAcceptor::initialize(prm_file);
1535 *  
1536 *   MinimalSurfaceProblem<2> minimal_surface_problem_2d;
1537 *   minimal_surface_problem_2d.run(parameters.formulation,
1538 *   parameters.tolerance);
1539 *   }
1540 *   catch (std::exception &exc)
1541 *   {
1542 *   std::cerr << std::endl
1543 *   << std::endl
1544 *   << "----------------------------------------------------"
1545 *   << std::endl;
1546 *   std::cerr << "Exception on processing: " << std::endl
1547 *   << exc.what() << std::endl
1548 *   << "Aborting!" << std::endl
1549 *   << "----------------------------------------------------"
1550 *   << std::endl;
1551 *  
1552 *   return 1;
1553 *   }
1554 *   catch (...)
1555 *   {
1556 *   std::cerr << std::endl
1557 *   << std::endl
1558 *   << "----------------------------------------------------"
1559 *   << std::endl;
1560 *   std::cerr << "Unknown exception!" << std::endl
1561 *   << "Aborting!" << std::endl
1562 *   << "----------------------------------------------------"
1563 *   << std::endl;
1564 *   return 1;
1565 *   }
1566 *   return 0;
1567 *   }
1568 * @endcode
1569<a name="step_72-Results"></a><h1>Results</h1>
1570
1571
1572Since there was no change to the physics of the problem that has first been analyzed
1573in @ref step_15 "step-15", there is nothing to report about that. The only outwardly noticeable
1574difference between them is that, by default, this program will only run 9 mesh
1575refinement steps (as opposed to @ref step_15 "step-15", which executes 11 refinements).
1576This will be observable in the simulation status that appears between the
1577header text that prints which assembly method is being used, and the final
1578timings. (All timings reported below were obtained in release mode.)
1579
1580@code
1581Mesh refinement step 0
1582 Initial residual: 1.53143
1583 Residual: 1.08746
1584 Residual: 0.966748
1585 Residual: 0.859602
1586 Residual: 0.766462
1587 Residual: 0.685475
1588
1589...
1590
1591Mesh refinement step 9
1592 Initial residual: 0.00924594
1593 Residual: 0.00831928
1594 Residual: 0.0074859
1595 Residual: 0.0067363
1596 Residual: 0.00606197
1597 Residual: 0.00545529
1598@endcode
1599
1600So what is interesting for us to compare is how long the assembly process takes
1601for the three different implementations, and to put that into some greater context.
1602Below is the output for the hand linearization (as computed on a circa 2012
1603four core, eight thread laptop -- but we're only really interested in the
1604relative time between the different implementations):
1605@code
1606******** Assembly approach ********
1607Unassisted implementation (full hand linearization).
1608
1609...
1610
1611+---------------------------------------------+------------+------------+
1612| Total wallclock time elapsed since start | 35.1s | |
1613| | | |
1614| Section | no. calls | wall time | % of total |
1615+---------------------------------+-----------+------------+------------+
1616| Assemble | 50 | 1.56s | 4.5% |
1617| Solve | 50 | 30.8s | 88% |
1618+---------------------------------+-----------+------------+------------+
1619@endcode
1620And for the implementation that linearizes the residual in an automated
1621manner using the Sacado dynamic forward AD number type:
1622@code
1623******** Assembly approach ********
1624Automated linearization of the finite element residual.
1625
1626...
1627
1628+---------------------------------------------+------------+------------+
1629| Total wallclock time elapsed since start | 40.1s | |
1630| | | |
1631| Section | no. calls | wall time | % of total |
1632+---------------------------------+-----------+------------+------------+
1633| Assemble | 50 | 8.8s | 22% |
1634| Solve | 50 | 28.6s | 71% |
1635+---------------------------------+-----------+------------+------------+
1636@endcode
1637And, lastly, for the implementation that computes both the residual and
1638its linearization directly from an energy functional (using nested Sacado
1639dynamic forward AD numbers):
1640@code
1641******** Assembly approach ********
1642Automated computation of finite element residual and linearization using a variational formulation.
1643
1644...
1645
1646+---------------------------------------------+------------+------------+
1647| Total wallclock time elapsed since start | 48.8s | |
1648| | | |
1649| Section | no. calls | wall time | % of total |
1650+---------------------------------+-----------+------------+------------+
1651| Assemble | 50 | 16.7s | 34% |
1652| Solve | 50 | 29.3s | 60% |
1653+---------------------------------+-----------+------------+------------+
1654@endcode
1655
1656It's evident that the more work that is passed off to the automatic differentiation
1657framework to perform, the more time is spent during the assembly process. Accumulated
1658over all refinement steps, using one level of automatic differentiation resulted
1659in @f$5.65 \times@f$ more computational time being spent in the assembly stage when
1660compared to unassisted assembly, while assembling the discrete linear system took
1661@f$10.7 \times@f$ longer when deriving directly from the energy functional.
1662Unsurprisingly, the overall time spent solving the linear system remained unchanged.
1663This means that the proportion of time spent in the solve phase to the assembly phase
1664shifted significantly as the number of times automated differentiation was performed
1665at the finite element level. For many, this might mean that leveraging higher-order
1666differentiation (at the finite element level) in production code leads to an
1667unacceptable overhead, but it may still be useful during the prototyping phase.
1668A good compromise between the two may, therefore, be the automated linearization
1669of the finite element residual, which offers a lot of convenience at a measurable,
1670but perhaps not unacceptable, cost. Alternatively, one could consider
1671not re-building the Newton matrix in every step -- a topic that is
1672explored in substantial depth in @ref step_77 "step-77".
1673
1674Of course, in practice the actual overhead is very much dependent on the problem being evaluated
1675(e.g., how many components there are in the solution, what the nature of the function
1676being differentiated is, etc.). So the exact results presented here should be
1677interpreted within the context of this scalar problem alone, and when it comes to
1678other problems, some preliminary investigation by the user is certainly warranted.
1679
1680
1681<a name="step_72-Possibilitiesforextensions"></a><h3> Possibilities for extensions </h3>
1682
1683
1684Like @ref step_71 "step-71", there are a few items related to automatic differentiation that could
1685be evaluated further:
1686- The use of other AD frameworks should be investigated, with the outlook that
1687 alternative implementations may provide performance benefits.
1688- It is also worth evaluating AD number types other than those that have been
1689 hard-coded into this tutorial. With regard to twice differentiable types
1690 employed at the finite-element level, mixed differentiation modes ("RAD-FAD")
1691 should in principle be more computationally efficient than the single
1692 mode ("FAD-FAD") types employed here. The reason that the RAD-FAD type was not
1693 selected by default is that, at the time of writing, there remain some
1694 bugs in its implementation within the Sacado library that lead to memory leaks.
1695 This is documented in the @ref auto_symb_diff topic.
1696- It might be the case that using reduced precision types (i.e., `float`) as the
1697 scalar types for the AD numbers could render a reduction in computational
1698 expense during assembly. Using `float` as the data type for the
1699 matrix and the residual is not unreasonable, given that the Newton
1700 update is only meant to get us closer to the solution, but not
1701 actually *to* the solution; as a consequence, it makes sense to
1702 consider using reduced-precision data types for computing these
1703 updates, and then accumulating these updates in a solution vector
1704 that uses the full `double` precision accuracy.
1705- One further method of possibly reducing resources during assembly is to frame
1706 the AD implementations as a constitutive model. This would be similar to the
1707 approach adopted in @ref step_71 "step-71", and pushes the starting point for the automatic
1708 differentiation one level higher up the chain of computations. This, in turn,
1709 means that less operations are tracked by the AD library, thereby reducing the
1710 cost of differentiating (even though one would perform the differentiation at
1711 each cell quadrature point).
1712- @ref step_77 "step-77" is yet another variation of @ref step_15 "step-15" that addresses a very
1713 different part of the problem: Line search and whether it is
1714 necessary to re-build the Newton matrix in every nonlinear
1715 iteration. Given that the results above show that using automatic
1716 differentiation comes at a cost, the techniques in @ref step_77 "step-77" have the
1717 potential to offset these costs to some degree. It would therefore
1718 be quite interesting to combine the current program with the
1719 techniques in @ref step_77 "step-77". For production codes, this would certainly be
1720 the way to go.
1721 *
1722 *
1723<a name="step_72-PlainProg"></a>
1724<h1> The plain program</h1>
1725@include "step-72.cc"
1726*/
void attach_dof_handler(const DoFHandler< dim, spacedim > &)
void distribute_dofs(const FiniteElement< dim, spacedim > &fe)
const unsigned int dofs_per_cell
Definition fe_data.h:436
static void estimate(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim - 1 > &quadrature, const std::map< types::boundary_id, const Function< spacedim, Number > * > &neumann_bc, const ReadVector< Number > &solution, Vector< float > &error, const ComponentMask &component_mask={}, const Function< spacedim > *coefficients=nullptr, const unsigned int n_threads=numbers::invalid_unsigned_int, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id, const types::material_id material_id=numbers::invalid_material_id, const Strategy strategy=cell_diameter_over_24)
static void initialize(const std::string &filename="", const std::string &output_filename="", const ParameterHandler::OutputStyle output_style_for_output_filename=ParameterHandler::Short, ParameterHandler &prm=ParameterAcceptor::prm, const ParameterHandler::OutputStyle output_style_for_filename=ParameterHandler::DefaultStyle)
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
@ wall_times
Definition timer.h:651
unsigned int n_active_cells() const
void refine_global(const unsigned int times=1)
virtual void execute_coarsening_and_refinement() override
Definition tria.cc:3320
virtual bool prepare_coarsening_and_refinement() override
Definition tria.cc:2805
float depth
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 reduction(Number *result, const Number *v, const size_type N)
__global__ void set(Number *val, const Number s, const size_type N)
#define AssertIndexRange(index, range)
#define AssertThrow(cond, exc)
void mesh_loop(const CellIteratorType &begin, const CellIteratorType &end, const CellWorkerFunctionType &cell_worker, const CopierType &copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const AssembleFlags flags=assemble_own_cells, const BoundaryWorkerFunctionType &boundary_worker=BoundaryWorkerFunctionType(), const FaceWorkerFunctionType &face_worker=FaceWorkerFunctionType(), const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition mesh_loop.h:281
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
@ update_JxW_values
Transformed quadrature weights.
@ update_gradients
Shape function gradients.
@ update_quadrature_points
Transformed quadrature points.
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
void hyper_ball(Triangulation< dim > &tria, const Point< dim > &center=Point< dim >(), const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
void refine_and_coarsen_fixed_number(Triangulation< dim, spacedim > &triangulation, const Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const unsigned int max_n_cells=std::numeric_limits< unsigned int >::max())
@ matrix
Contents is actually a matrix.
void cell_matrix(FullMatrix< double > &M, const FEValuesBase< dim > &fe, const FEValuesBase< dim > &fetest, const ArrayView< const std::vector< double > > &velocity, const double factor=1.)
Definition advection.h:74
void cell_residual(Vector< double > &result, const FEValuesBase< dim > &fe, const std::vector< Tensor< 1, dim > > &input, const ArrayView< const std::vector< double > > &velocity, double factor=1.)
Definition advection.h:130
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition utilities.cc:191
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
constexpr const ReferenceCell Line
constexpr ReturnType< rank, T >::value_type & extract(T &t, const ArrayType &indices)
std::string int_to_string(const unsigned int value, const unsigned int digits=numbers::invalid_unsigned_int)
Definition utilities.cc:470
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})
void 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)
void copy(const T *begin, const T *end, U *dest)
int(& functions)(const void *v1, const void *v2)
static constexpr double PI
Definition numbers.h:259
STL namespace.
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition types.h:32
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation