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-55.h
Go to the documentation of this file.
1 
234  *
235  * namespace LA
236  * {
237  * #if defined(DEAL_II_WITH_PETSC) && !defined(DEAL_II_PETSC_WITH_COMPLEX) && \
238  * !(defined(DEAL_II_WITH_TRILINOS) && defined(FORCE_USE_OF_TRILINOS))
239  * using namespace dealii::LinearAlgebraPETSc;
240  * # define USE_PETSC_LA
241  * #elif defined(DEAL_II_WITH_TRILINOS)
242  * using namespace dealii::LinearAlgebraTrilinos;
243  * #else
244  * # error DEAL_II_WITH_PETSC or DEAL_II_WITH_TRILINOS required
245  * #endif
246  * } // namespace LA
247  *
248  * #include <deal.II/lac/vector.h>
249  * #include <deal.II/lac/full_matrix.h>
250  * #include <deal.II/lac/solver_cg.h>
251  * #include <deal.II/lac/solver_gmres.h>
252  * #include <deal.II/lac/solver_minres.h>
253  * #include <deal.II/lac/affine_constraints.h>
254  * #include <deal.II/lac/dynamic_sparsity_pattern.h>
255  *
256  * #include <deal.II/lac/petsc_sparse_matrix.h>
257  * #include <deal.II/lac/petsc_vector.h>
258  * #include <deal.II/lac/petsc_solver.h>
259  * #include <deal.II/lac/petsc_precondition.h>
260  *
261  * #include <deal.II/grid/grid_generator.h>
262  * #include <deal.II/grid/tria_accessor.h>
263  * #include <deal.II/grid/tria_iterator.h>
264  * #include <deal.II/grid/manifold_lib.h>
265  * #include <deal.II/grid/grid_tools.h>
266  * #include <deal.II/dofs/dof_handler.h>
267  * #include <deal.II/dofs/dof_renumbering.h>
268  * #include <deal.II/dofs/dof_accessor.h>
269  * #include <deal.II/dofs/dof_tools.h>
270  * #include <deal.II/fe/fe_values.h>
271  * #include <deal.II/fe/fe_q.h>
272  * #include <deal.II/fe/fe_system.h>
273  * #include <deal.II/numerics/vector_tools.h>
274  * #include <deal.II/numerics/data_out.h>
275  * #include <deal.II/numerics/error_estimator.h>
276  *
277  * #include <deal.II/base/utilities.h>
278  * #include <deal.II/base/conditional_ostream.h>
279  * #include <deal.II/base/index_set.h>
280  * #include <deal.II/lac/sparsity_tools.h>
281  * #include <deal.II/distributed/tria.h>
282  * #include <deal.II/distributed/grid_refinement.h>
283  *
284  * #include <cmath>
285  * #include <fstream>
286  * #include <iostream>
287  *
288  * namespace Step55
289  * {
290  * using namespace dealii;
291  *
292  * @endcode
293  *
294  *
295  * <a name="Linearsolversandpreconditioners"></a>
296  * <h3>Linear solvers and preconditioners</h3>
297  *
298 
299  *
300  * We need a few helper classes to represent our solver strategy
301  * described in the introduction.
302  *
303 
304  *
305  *
306  * @code
307  * namespace LinearSolvers
308  * {
309  * @endcode
310  *
311  * This class exposes the action of applying the inverse of a giving
312  * matrix via the function InverseMatrix::vmult(). Internally, the
313  * inverse is not formed explicitly. Instead, a linear solver with CG
314  * is performed. This class extends the InverseMatrix class in @ref step_22 "step-22"
315  * with an option to specify a preconditioner, and to allow for different
316  * vector types in the vmult function.
317  *
318  * @code
319  * template <class Matrix, class Preconditioner>
320  * class InverseMatrix : public Subscriptor
321  * {
322  * public:
323  * InverseMatrix(const Matrix &m, const Preconditioner &preconditioner);
324  *
325  * template <typename VectorType>
326  * void vmult(VectorType &dst, const VectorType &src) const;
327  *
328  * private:
330  * const Preconditioner & preconditioner;
331  * };
332  *
333  *
334  * template <class Matrix, class Preconditioner>
335  * InverseMatrix<Matrix, Preconditioner>::InverseMatrix(
336  * const Matrix & m,
337  * const Preconditioner &preconditioner)
338  * : matrix(&m)
339  * , preconditioner(preconditioner)
340  * {}
341  *
342  *
343  *
344  * template <class Matrix, class Preconditioner>
345  * template <typename VectorType>
346  * void
347  * InverseMatrix<Matrix, Preconditioner>::vmult(VectorType & dst,
348  * const VectorType &src) const
349  * {
350  * SolverControl solver_control(src.size(), 1e-8 * src.l2_norm());
351  * SolverCG<LA::MPI::Vector> cg(solver_control);
352  * dst = 0;
353  *
354  * try
355  * {
356  * cg.solve(*matrix, dst, src, preconditioner);
357  * }
358  * catch (std::exception &e)
359  * {
360  * Assert(false, ExcMessage(e.what()));
361  * }
362  * }
363  *
364  *
365  * @endcode
366  *
367  * The class A template class for a simple block diagonal preconditioner
368  * for 2x2 matrices.
369  *
370  * @code
371  * template <class PreconditionerA, class PreconditionerS>
372  * class BlockDiagonalPreconditioner : public Subscriptor
373  * {
374  * public:
375  * BlockDiagonalPreconditioner(const PreconditionerA &preconditioner_A,
376  * const PreconditionerS &preconditioner_S);
377  *
378  * void vmult(LA::MPI::BlockVector & dst,
379  * const LA::MPI::BlockVector &src) const;
380  *
381  * private:
382  * const PreconditionerA &preconditioner_A;
383  * const PreconditionerS &preconditioner_S;
384  * };
385  *
386  * template <class PreconditionerA, class PreconditionerS>
387  * BlockDiagonalPreconditioner<PreconditionerA, PreconditionerS>::
388  * BlockDiagonalPreconditioner(const PreconditionerA &preconditioner_A,
389  * const PreconditionerS &preconditioner_S)
390  * : preconditioner_A(preconditioner_A)
391  * , preconditioner_S(preconditioner_S)
392  * {}
393  *
394  *
395  * template <class PreconditionerA, class PreconditionerS>
396  * void BlockDiagonalPreconditioner<PreconditionerA, PreconditionerS>::vmult(
397  * LA::MPI::BlockVector & dst,
398  * const LA::MPI::BlockVector &src) const
399  * {
400  * preconditioner_A.vmult(dst.block(0), src.block(0));
401  * preconditioner_S.vmult(dst.block(1), src.block(1));
402  * }
403  *
404  * } // namespace LinearSolvers
405  *
406  * @endcode
407  *
408  *
409  * <a name="Problemsetup"></a>
410  * <h3>Problem setup</h3>
411  *
412 
413  *
414  * The following classes represent the right hand side and the exact
415  * solution for the test problem.
416  *
417 
418  *
419  *
420  * @code
421  * template <int dim>
422  * class RightHandSide : public Function<dim>
423  * {
424  * public:
425  * RightHandSide()
426  * : Function<dim>(dim + 1)
427  * {}
428  *
429  * virtual void vector_value(const Point<dim> &p,
430  * Vector<double> & value) const override;
431  * };
432  *
433  *
434  * template <int dim>
435  * void RightHandSide<dim>::vector_value(const Point<dim> &p,
436  * Vector<double> & values) const
437  * {
438  * const double R_x = p[0];
439  * const double R_y = p[1];
440  *
441  * const double pi = numbers::PI;
442  * const double pi2 = pi * pi;
443  * values[0] =
444  * -1.0L / 2.0L * (-2 * sqrt(25.0 + 4 * pi2) + 10.0) *
445  * exp(R_x * (-2 * sqrt(25.0 + 4 * pi2) + 10.0)) -
446  * 0.4 * pi2 * exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * cos(2 * R_y * pi) +
447  * 0.1 * pow(-sqrt(25.0 + 4 * pi2) + 5.0, 2) *
448  * exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * cos(2 * R_y * pi);
449  * values[1] = 0.2 * pi * (-sqrt(25.0 + 4 * pi2) + 5.0) *
450  * exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * sin(2 * R_y * pi) -
451  * 0.05 * pow(-sqrt(25.0 + 4 * pi2) + 5.0, 3) *
452  * exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * sin(2 * R_y * pi) /
453  * pi;
454  * values[2] = 0;
455  * }
456  *
457  *
458  * template <int dim>
459  * class ExactSolution : public Function<dim>
460  * {
461  * public:
462  * ExactSolution()
463  * : Function<dim>(dim + 1)
464  * {}
465  *
466  * virtual void vector_value(const Point<dim> &p,
467  * Vector<double> & value) const override;
468  * };
469  *
470  * template <int dim>
471  * void ExactSolution<dim>::vector_value(const Point<dim> &p,
472  * Vector<double> & values) const
473  * {
474  * const double R_x = p[0];
475  * const double R_y = p[1];
476  *
477  * const double pi = numbers::PI;
478  * const double pi2 = pi * pi;
479  * values[0] =
480  * -exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * cos(2 * R_y * pi) + 1;
481  * values[1] = (1.0L / 2.0L) * (-sqrt(25.0 + 4 * pi2) + 5.0) *
482  * exp(R_x * (-sqrt(25.0 + 4 * pi2) + 5.0)) * sin(2 * R_y * pi) /
483  * pi;
484  * values[2] =
485  * -1.0L / 2.0L * exp(R_x * (-2 * sqrt(25.0 + 4 * pi2) + 10.0)) -
486  * 2.0 *
487  * (-6538034.74494422 +
488  * 0.0134758939981709 * exp(4 * sqrt(25.0 + 4 * pi2))) /
489  * (-80.0 * exp(3 * sqrt(25.0 + 4 * pi2)) +
490  * 16.0 * sqrt(25.0 + 4 * pi2) * exp(3 * sqrt(25.0 + 4 * pi2))) -
491  * 1634508.68623606 * exp(-3.0 * sqrt(25.0 + 4 * pi2)) /
492  * (-10.0 + 2.0 * sqrt(25.0 + 4 * pi2)) +
493  * (-0.00673794699908547 * exp(sqrt(25.0 + 4 * pi2)) +
494  * 3269017.37247211 * exp(-3 * sqrt(25.0 + 4 * pi2))) /
495  * (-8 * sqrt(25.0 + 4 * pi2) + 40.0) +
496  * 0.00336897349954273 * exp(1.0 * sqrt(25.0 + 4 * pi2)) /
497  * (-10.0 + 2.0 * sqrt(25.0 + 4 * pi2));
498  * }
499  *
500  *
501  *
502  * @endcode
503  *
504  *
505  * <a name="Themainprogram"></a>
506  * <h3>The main program</h3>
507  *
508 
509  *
510  * The main class is very similar to @ref step_40 "step-40", except that matrices and
511  * vectors are now block versions, and we store a std::vector<IndexSet>
512  * for owned and relevant DoFs instead of a single IndexSet. We have
513  * exactly two IndexSets, one for all velocity unknowns and one for all
514  * pressure unknowns.
515  *
516  * @code
517  * template <int dim>
518  * class StokesProblem
519  * {
520  * public:
521  * StokesProblem(unsigned int velocity_degree);
522  *
523  * void run();
524  *
525  * private:
526  * void make_grid();
527  * void setup_system();
528  * void assemble_system();
529  * void solve();
530  * void refine_grid();
531  * void output_results(const unsigned int cycle) const;
532  *
533  * unsigned int velocity_degree;
534  * double viscosity;
535  * MPI_Comm mpi_communicator;
536  *
537  * FESystem<dim> fe;
539  * DoFHandler<dim> dof_handler;
540  *
541  * std::vector<IndexSet> owned_partitioning;
542  * std::vector<IndexSet> relevant_partitioning;
543  *
544  * AffineConstraints<double> constraints;
545  *
546  * LA::MPI::BlockSparseMatrix system_matrix;
547  * LA::MPI::BlockSparseMatrix preconditioner_matrix;
548  * LA::MPI::BlockVector locally_relevant_solution;
549  * LA::MPI::BlockVector system_rhs;
550  *
551  * ConditionalOStream pcout;
552  * TimerOutput computing_timer;
553  * };
554  *
555  *
556  *
557  * template <int dim>
558  * StokesProblem<dim>::StokesProblem(unsigned int velocity_degree)
559  * : velocity_degree(velocity_degree)
560  * , viscosity(0.1)
561  * , mpi_communicator(MPI_COMM_WORLD)
562  * , fe(FE_Q<dim>(velocity_degree), dim, FE_Q<dim>(velocity_degree - 1), 1)
563  * , triangulation(mpi_communicator,
567  * , dof_handler(triangulation)
568  * , pcout(std::cout,
569  * (Utilities::MPI::this_mpi_process(mpi_communicator) == 0))
570  * , computing_timer(mpi_communicator,
571  * pcout,
574  * {}
575  *
576  *
577  * @endcode
578  *
579  * The Kovasnay flow is defined on the domain [-0.5, 1.5]^2, which we
580  * create by passing the min and max values to GridGenerator::hyper_cube.
581  *
582  * @code
583  * template <int dim>
584  * void StokesProblem<dim>::make_grid()
585  * {
587  * triangulation.refine_global(3);
588  * }
589  *
590  * @endcode
591  *
592  *
593  * <a name="SystemSetup"></a>
594  * <h3>System Setup</h3>
595  *
596 
597  *
598  * The construction of the block matrices and vectors is new compared to
599  * @ref step_40 "step-40" and is different compared to serial codes like @ref step_22 "step-22", because
600  * we need to supply the set of rows that belong to our processor.
601  *
602  * @code
603  * template <int dim>
604  * void StokesProblem<dim>::setup_system()
605  * {
606  * TimerOutput::Scope t(computing_timer, "setup");
607  *
608  * dof_handler.distribute_dofs(fe);
609  *
610  * @endcode
611  *
612  * Put all dim velocities into block 0 and the pressure into block 1,
613  * then reorder the unknowns by block. Finally count how many unknowns
614  * we have per block.
615  *
616  * @code
617  * std::vector<unsigned int> stokes_sub_blocks(dim + 1, 0);
618  * stokes_sub_blocks[dim] = 1;
619  * DoFRenumbering::component_wise(dof_handler, stokes_sub_blocks);
620  *
621  * const std::vector<types::global_dof_index> dofs_per_block =
622  * DoFTools::count_dofs_per_fe_block(dof_handler, stokes_sub_blocks);
623  *
624  * const unsigned int n_u = dofs_per_block[0];
625  * const unsigned int n_p = dofs_per_block[1];
626  *
627  * pcout << " Number of degrees of freedom: " << dof_handler.n_dofs() << " ("
628  * << n_u << '+' << n_p << ')' << std::endl;
629  *
630  * @endcode
631  *
632  * We split up the IndexSet for locally owned and locally relevant DoFs
633  * into two IndexSets based on how we want to create the block matrices
634  * and vectors.
635  *
636  * @code
637  * owned_partitioning.resize(2);
638  * owned_partitioning[0] = dof_handler.locally_owned_dofs().get_view(0, n_u);
639  * owned_partitioning[1] =
640  * dof_handler.locally_owned_dofs().get_view(n_u, n_u + n_p);
641  *
642  * IndexSet locally_relevant_dofs;
643  * DoFTools::extract_locally_relevant_dofs(dof_handler, locally_relevant_dofs);
644  * relevant_partitioning.resize(2);
645  * relevant_partitioning[0] = locally_relevant_dofs.get_view(0, n_u);
646  * relevant_partitioning[1] = locally_relevant_dofs.get_view(n_u, n_u + n_p);
647  *
648  * @endcode
649  *
650  * Setting up the constraints for boundary conditions and hanging nodes
651  * is identical to @ref step_40 "step-40". Rven though we don't have any hanging nodes
652  * because we only perform global refinement, it is still a good idea
653  * to put this function call in, in case adaptive refinement gets
654  * introduced later.
655  *
656  * @code
657  * {
658  * constraints.reinit(locally_relevant_dofs);
659  *
660  * FEValuesExtractors::Vector velocities(0);
661  * DoFTools::make_hanging_node_constraints(dof_handler, constraints);
662  * VectorTools::interpolate_boundary_values(dof_handler,
663  * 0,
664  * ExactSolution<dim>(),
665  * constraints,
666  * fe.component_mask(velocities));
667  * constraints.close();
668  * }
669  *
670  * @endcode
671  *
672  * Now we create the system matrix based on a BlockDynamicSparsityPattern.
673  * We know that we won't have coupling between different velocity
674  * components (because we use the laplace and not the deformation tensor)
675  * and no coupling between pressure with its test functions, so we use
676  * a Table to communicate this coupling information to
678  *
679  * @code
680  * {
681  * system_matrix.clear();
682  *
683  * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
684  * for (unsigned int c = 0; c < dim + 1; ++c)
685  * for (unsigned int d = 0; d < dim + 1; ++d)
686  * if (c == dim && d == dim)
687  * coupling[c][d] = DoFTools::none;
688  * else if (c == dim || d == dim || c == d)
689  * coupling[c][d] = DoFTools::always;
690  * else
691  * coupling[c][d] = DoFTools::none;
692  *
693  * BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
694  *
696  * dof_handler, coupling, dsp, constraints, false);
697  *
699  * dsp,
700  * dof_handler.locally_owned_dofs(),
701  * mpi_communicator,
702  * locally_relevant_dofs);
703  *
704  * system_matrix.reinit(owned_partitioning, dsp, mpi_communicator);
705  * }
706  *
707  * @endcode
708  *
709  * The preconditioner matrix has a different coupling (we only fill in
710  * the 1,1 block with the mass matrix), otherwise this code is identical
711  * to the construction of the system_matrix above.
712  *
713  * @code
714  * {
715  * preconditioner_matrix.clear();
716  *
717  * Table<2, DoFTools::Coupling> coupling(dim + 1, dim + 1);
718  * for (unsigned int c = 0; c < dim + 1; ++c)
719  * for (unsigned int d = 0; d < dim + 1; ++d)
720  * if (c == dim && d == dim)
721  * coupling[c][d] = DoFTools::always;
722  * else
723  * coupling[c][d] = DoFTools::none;
724  *
725  * BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
726  *
728  * dof_handler, coupling, dsp, constraints, false);
730  * dsp,
731  * Utilities::MPI::all_gather(mpi_communicator,
732  * dof_handler.locally_owned_dofs()),
733  * mpi_communicator,
734  * locally_relevant_dofs);
735  * preconditioner_matrix.reinit(owned_partitioning,
736  * @endcode
737  *
738  * owned_partitioning,
739  *
740  * @code
741  * dsp,
742  * mpi_communicator);
743  * }
744  *
745  * @endcode
746  *
747  * Finally, we construct the block vectors with the right sizes. The
748  * function call with two std::vector<IndexSet> will create a ghosted
749  * vector.
750  *
751  * @code
752  * locally_relevant_solution.reinit(owned_partitioning,
753  * relevant_partitioning,
754  * mpi_communicator);
755  * system_rhs.reinit(owned_partitioning, mpi_communicator);
756  * }
757  *
758  *
759  *
760  * @endcode
761  *
762  *
763  * <a name="Assembly"></a>
764  * <h3>Assembly</h3>
765  *
766 
767  *
768  * This function assembles the system matrix, the preconditioner matrix,
769  * and the right hand side. The code is pretty standard.
770  *
771  * @code
772  * template <int dim>
773  * void StokesProblem<dim>::assemble_system()
774  * {
775  * TimerOutput::Scope t(computing_timer, "assembly");
776  *
777  * system_matrix = 0;
778  * preconditioner_matrix = 0;
779  * system_rhs = 0;
780  *
781  * const QGauss<dim> quadrature_formula(velocity_degree + 1);
782  *
783  * FEValues<dim> fe_values(fe,
784  * quadrature_formula,
787  *
788  * const unsigned int dofs_per_cell = fe.dofs_per_cell;
789  * const unsigned int n_q_points = quadrature_formula.size();
790  *
791  * FullMatrix<double> cell_matrix(dofs_per_cell, dofs_per_cell);
792  * FullMatrix<double> cell_matrix2(dofs_per_cell, dofs_per_cell);
793  * Vector<double> cell_rhs(dofs_per_cell);
794  *
795  * const RightHandSide<dim> right_hand_side;
796  * std::vector<Vector<double>> rhs_values(n_q_points, Vector<double>(dim + 1));
797  *
798  * std::vector<Tensor<2, dim>> grad_phi_u(dofs_per_cell);
799  * std::vector<double> div_phi_u(dofs_per_cell);
800  * std::vector<double> phi_p(dofs_per_cell);
801  *
802  * std::vector<types::global_dof_index> local_dof_indices(dofs_per_cell);
803  * const FEValuesExtractors::Vector velocities(0);
804  * const FEValuesExtractors::Scalar pressure(dim);
805  *
806  * for (const auto &cell : dof_handler.active_cell_iterators())
807  * if (cell->is_locally_owned())
808  * {
809  * cell_matrix = 0;
810  * cell_matrix2 = 0;
811  * cell_rhs = 0;
812  *
813  * fe_values.reinit(cell);
814  * right_hand_side.vector_value_list(fe_values.get_quadrature_points(),
815  * rhs_values);
816  * for (unsigned int q = 0; q < n_q_points; ++q)
817  * {
818  * for (unsigned int k = 0; k < dofs_per_cell; ++k)
819  * {
820  * grad_phi_u[k] = fe_values[velocities].gradient(k, q);
821  * div_phi_u[k] = fe_values[velocities].divergence(k, q);
822  * phi_p[k] = fe_values[pressure].value(k, q);
823  * }
824  *
825  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
826  * {
827  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
828  * {
829  * cell_matrix(i, j) +=
830  * (viscosity *
831  * scalar_product(grad_phi_u[i], grad_phi_u[j]) -
832  * div_phi_u[i] * phi_p[j] - phi_p[i] * div_phi_u[j]) *
833  * fe_values.JxW(q);
834  *
835  * cell_matrix2(i, j) += 1.0 / viscosity * phi_p[i] *
836  * phi_p[j] * fe_values.JxW(q);
837  * }
838  *
839  * const unsigned int component_i =
840  * fe.system_to_component_index(i).first;
841  * cell_rhs(i) += fe_values.shape_value(i, q) *
842  * rhs_values[q](component_i) * fe_values.JxW(q);
843  * }
844  * }
845  *
846  *
847  * cell->get_dof_indices(local_dof_indices);
848  * constraints.distribute_local_to_global(cell_matrix,
849  * cell_rhs,
850  * local_dof_indices,
851  * system_matrix,
852  * system_rhs);
853  *
854  * constraints.distribute_local_to_global(cell_matrix2,
855  * local_dof_indices,
856  * preconditioner_matrix);
857  * }
858  *
859  * system_matrix.compress(VectorOperation::add);
860  * preconditioner_matrix.compress(VectorOperation::add);
861  * system_rhs.compress(VectorOperation::add);
862  * }
863  *
864  *
865  *
866  * @endcode
867  *
868  *
869  * <a name="Solving"></a>
870  * <h3>Solving</h3>
871  *
872 
873  *
874  * This function solves the linear system with MINRES with a block diagonal
875  * preconditioner and AMG for the two diagonal blocks as described in the
876  * introduction. The preconditioner applies a v cycle to the 0,0 block
877  * and a CG with the mass matrix for the 1,1 block (the Schur complement).
878  *
879  * @code
880  * template <int dim>
881  * void StokesProblem<dim>::solve()
882  * {
883  * TimerOutput::Scope t(computing_timer, "solve");
884  *
885  * LA::MPI::PreconditionAMG prec_A;
886  * {
887  * LA::MPI::PreconditionAMG::AdditionalData data;
888  *
889  * #ifdef USE_PETSC_LA
890  * data.symmetric_operator = true;
891  * #endif
892  * prec_A.initialize(system_matrix.block(0, 0), data);
893  * }
894  *
895  * LA::MPI::PreconditionAMG prec_S;
896  * {
897  * LA::MPI::PreconditionAMG::AdditionalData data;
898  *
899  * #ifdef USE_PETSC_LA
900  * data.symmetric_operator = true;
901  * #endif
902  * prec_S.initialize(preconditioner_matrix.block(1, 1), data);
903  * }
904  *
905  * @endcode
906  *
907  * The InverseMatrix is used to solve for the mass matrix:
908  *
909  * @code
910  * using mp_inverse_t = LinearSolvers::InverseMatrix<LA::MPI::SparseMatrix,
912  * const mp_inverse_t mp_inverse(preconditioner_matrix.block(1, 1), prec_S);
913  *
914  * @endcode
915  *
916  * This constructs the block preconditioner based on the preconditioners
917  * for the individual blocks defined above.
918  *
919  * @code
920  * const LinearSolvers::BlockDiagonalPreconditioner<LA::MPI::PreconditionAMG,
921  * mp_inverse_t>
922  * preconditioner(prec_A, mp_inverse);
923  *
924  * @endcode
925  *
926  * With that, we can finally set up the linear solver and solve the system:
927  *
928  * @code
929  * SolverControl solver_control(system_matrix.m(),
930  * 1e-10 * system_rhs.l2_norm());
931  *
932  * SolverMinRes<LA::MPI::BlockVector> solver(solver_control);
933  *
934  * LA::MPI::BlockVector distributed_solution(owned_partitioning,
935  * mpi_communicator);
936  *
937  * constraints.set_zero(distributed_solution);
938  *
939  * solver.solve(system_matrix,
940  * distributed_solution,
941  * system_rhs,
942  * preconditioner);
943  *
944  * pcout << " Solved in " << solver_control.last_step() << " iterations."
945  * << std::endl;
946  *
947  * constraints.distribute(distributed_solution);
948  *
949  * @endcode
950  *
951  * Like in @ref step_56 "step-56", we subtract the mean pressure to allow error
952  * computations against our reference solution, which has a mean value
953  * of zero.
954  *
955  * @code
956  * locally_relevant_solution = distributed_solution;
957  * const double mean_pressure =
958  * VectorTools::compute_mean_value(dof_handler,
959  * QGauss<dim>(velocity_degree + 2),
960  * locally_relevant_solution,
961  * dim);
962  * distributed_solution.block(1).add(-mean_pressure);
963  * locally_relevant_solution.block(1) = distributed_solution.block(1);
964  * }
965  *
966  *
967  *
968  * @endcode
969  *
970  *
971  * <a name="Therest"></a>
972  * <h3>The rest</h3>
973  *
974 
975  *
976  * The remainder of the code that deals with mesh refinement, output, and
977  * the main loop is pretty standard.
978  *
979  * @code
980  * template <int dim>
981  * void StokesProblem<dim>::refine_grid()
982  * {
983  * TimerOutput::Scope t(computing_timer, "refine");
984  *
985  * triangulation.refine_global();
986  * }
987  *
988  *
989  *
990  * template <int dim>
991  * void StokesProblem<dim>::output_results(const unsigned int cycle) const
992  * {
993  * {
994  * const ComponentSelectFunction<dim> pressure_mask(dim, dim + 1);
995  * const ComponentSelectFunction<dim> velocity_mask(std::make_pair(0, dim),
996  * dim + 1);
997  *
998  * Vector<double> cellwise_errors(triangulation.n_active_cells());
999  * QGauss<dim> quadrature(velocity_degree + 2);
1000  *
1001  * VectorTools::integrate_difference(dof_handler,
1002  * locally_relevant_solution,
1003  * ExactSolution<dim>(),
1004  * cellwise_errors,
1005  * quadrature,
1007  * &velocity_mask);
1008  *
1009  * const double error_u_l2 =
1011  * cellwise_errors,
1013  *
1014  * VectorTools::integrate_difference(dof_handler,
1015  * locally_relevant_solution,
1016  * ExactSolution<dim>(),
1017  * cellwise_errors,
1018  * quadrature,
1020  * &pressure_mask);
1021  *
1022  * const double error_p_l2 =
1024  * cellwise_errors,
1026  *
1027  * pcout << "error: u_0: " << error_u_l2 << " p_0: " << error_p_l2
1028  * << std::endl;
1029  * }
1030  *
1031  *
1032  * std::vector<std::string> solution_names(dim, "velocity");
1033  * solution_names.emplace_back("pressure");
1034  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
1035  * data_component_interpretation(
1037  * data_component_interpretation.push_back(
1039  *
1040  * DataOut<dim> data_out;
1041  * data_out.attach_dof_handler(dof_handler);
1042  * data_out.add_data_vector(locally_relevant_solution,
1043  * solution_names,
1045  * data_component_interpretation);
1046  *
1047  * LA::MPI::BlockVector interpolated;
1048  * interpolated.reinit(owned_partitioning, MPI_COMM_WORLD);
1049  * VectorTools::interpolate(dof_handler, ExactSolution<dim>(), interpolated);
1050  *
1051  * LA::MPI::BlockVector interpolated_relevant(owned_partitioning,
1052  * relevant_partitioning,
1053  * MPI_COMM_WORLD);
1054  * interpolated_relevant = interpolated;
1055  * {
1056  * std::vector<std::string> solution_names(dim, "ref_u");
1057  * solution_names.emplace_back("ref_p");
1058  * data_out.add_data_vector(interpolated_relevant,
1059  * solution_names,
1061  * data_component_interpretation);
1062  * }
1063  *
1064  *
1065  * Vector<float> subdomain(triangulation.n_active_cells());
1066  * for (unsigned int i = 0; i < subdomain.size(); ++i)
1067  * subdomain(i) = triangulation.locally_owned_subdomain();
1068  * data_out.add_data_vector(subdomain, "subdomain");
1069  *
1070  * data_out.build_patches();
1071  *
1072  * data_out.write_vtu_with_pvtu_record(
1073  * "./", "solution", cycle, mpi_communicator, 2);
1074  * }
1075  *
1076  *
1077  *
1078  * template <int dim>
1079  * void StokesProblem<dim>::run()
1080  * {
1081  * #ifdef USE_PETSC_LA
1082  * pcout << "Running using PETSc." << std::endl;
1083  * #else
1084  * pcout << "Running using Trilinos." << std::endl;
1085  * #endif
1086  * const unsigned int n_cycles = 5;
1087  * for (unsigned int cycle = 0; cycle < n_cycles; ++cycle)
1088  * {
1089  * pcout << "Cycle " << cycle << ':' << std::endl;
1090  *
1091  * if (cycle == 0)
1092  * make_grid();
1093  * else
1094  * refine_grid();
1095  *
1096  * setup_system();
1097  *
1098  * assemble_system();
1099  * solve();
1100  *
1101  * if (Utilities::MPI::n_mpi_processes(mpi_communicator) <= 32)
1102  * {
1103  * TimerOutput::Scope t(computing_timer, "output");
1104  * output_results(cycle);
1105  * }
1106  *
1107  * computing_timer.print_summary();
1108  * computing_timer.reset();
1109  *
1110  * pcout << std::endl;
1111  * }
1112  * }
1113  * } // namespace Step55
1114  *
1115  *
1116  *
1117  * int main(int argc, char *argv[])
1118  * {
1119  * try
1120  * {
1121  * using namespace dealii;
1122  * using namespace Step55;
1123  *
1124  * Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, 1);
1125  *
1126  * StokesProblem<2> problem(2);
1127  * problem.run();
1128  * }
1129  * catch (std::exception &exc)
1130  * {
1131  * std::cerr << std::endl
1132  * << std::endl
1133  * << "----------------------------------------------------"
1134  * << std::endl;
1135  * std::cerr << "Exception on processing: " << std::endl
1136  * << exc.what() << std::endl
1137  * << "Aborting!" << std::endl
1138  * << "----------------------------------------------------"
1139  * << std::endl;
1140  *
1141  * return 1;
1142  * }
1143  * catch (...)
1144  * {
1145  * std::cerr << std::endl
1146  * << std::endl
1147  * << "----------------------------------------------------"
1148  * << std::endl;
1149  * std::cerr << "Unknown exception!" << std::endl
1150  * << "Aborting!" << std::endl
1151  * << "----------------------------------------------------"
1152  * << std::endl;
1153  * return 1;
1154  * }
1155  *
1156  * return 0;
1157  * }
1158  * @endcode
1159 <a name="Results"></a><h1>Results</h1>
1160 
1161 
1162 As expected from the discussion above, the number of iterations is independent
1163 of the number of processors and only very slightly dependent on @f$h@f$:
1164 
1165 <table>
1166 <tr>
1167  <th colspan="2">PETSc</th>
1168  <th colspan="8">number of processors</th>
1169 </tr>
1170 <tr>
1171  <th>cycle</th>
1172  <th>dofs</th>
1173  <th>1</th>
1174  <th>2</th>
1175  <th>4</th>
1176  <th>8</th>
1177  <th>16</th>
1178  <th>32</th>
1179  <th>64</th>
1180  <th>128</th>
1181 </tr>
1182 <tr>
1183  <td>0</td>
1184  <td>659</td>
1185  <td>49</td>
1186  <td>49</td>
1187  <td>49</td>
1188  <td>51</td>
1189  <td>51</td>
1190  <td>51</td>
1191  <td>49</td>
1192  <td>49</td>
1193 </tr>
1194 <tr>
1195  <td>1</td>
1196  <td>2467</td>
1197  <td>52</td>
1198  <td>52</td>
1199  <td>52</td>
1200  <td>52</td>
1201  <td>52</td>
1202  <td>54</td>
1203  <td>54</td>
1204  <td>53</td>
1205 </tr>
1206 <tr>
1207  <td>2</td>
1208  <td>9539</td>
1209  <td>56</td>
1210  <td>56</td>
1211  <td>56</td>
1212  <td>54</td>
1213  <td>56</td>
1214  <td>56</td>
1215  <td>54</td>
1216  <td>56</td>
1217 </tr>
1218 <tr>
1219  <td>3</td>
1220  <td>37507</td>
1221  <td>57</td>
1222  <td>57</td>
1223  <td>57</td>
1224  <td>57</td>
1225  <td>57</td>
1226  <td>56</td>
1227  <td>57</td>
1228  <td>56</td>
1229 </tr>
1230 <tr>
1231  <td>4</td>
1232  <td>148739</td>
1233  <td>58</td>
1234  <td>59</td>
1235  <td>57</td>
1236  <td>59</td>
1237  <td>57</td>
1238  <td>57</td>
1239  <td>57</td>
1240  <td>57</td>
1241 </tr>
1242 <tr>
1243  <td>5</td>
1244  <td>592387</td>
1245  <td>60</td>
1246  <td>60</td>
1247  <td>59</td>
1248  <td>59</td>
1249  <td>59</td>
1250  <td>59</td>
1251  <td>59</td>
1252  <td>59</td>
1253 </tr>
1254 <tr>
1255  <td>6</td>
1256  <td>2364419</td>
1257  <td>62</td>
1258  <td>62</td>
1259  <td>61</td>
1260  <td>61</td>
1261  <td>61</td>
1262  <td>61</td>
1263  <td>61</td>
1264  <td>61</td>
1265 </tr>
1266 </table>
1267 
1268 <table>
1269 <tr>
1270  <th colspan="2">Trilinos</th>
1271  <th colspan="8">number of processors</th>
1272 </tr>
1273 <tr>
1274  <th>cycle</th>
1275  <th>dofs</th>
1276  <th>1</th>
1277  <th>2</th>
1278  <th>4</th>
1279  <th>8</th>
1280  <th>16</th>
1281  <th>32</th>
1282  <th>64</th>
1283  <th>128</th>
1284 </tr>
1285 <tr>
1286  <td>0</td>
1287  <td>659</td>
1288  <td>37</td>
1289  <td>37</td>
1290  <td>37</td>
1291  <td>37</td>
1292  <td>37</td>
1293  <td>37</td>
1294  <td>37</td>
1295  <td>37</td>
1296 </tr>
1297 <tr>
1298  <td>1</td>
1299  <td>2467</td>
1300  <td>92</td>
1301  <td>89</td>
1302  <td>89</td>
1303  <td>82</td>
1304  <td>86</td>
1305  <td>81</td>
1306  <td>78</td>
1307  <td>78</td>
1308 </tr>
1309 <tr>
1310  <td>2</td>
1311  <td>9539</td>
1312  <td>102</td>
1313  <td>99</td>
1314  <td>96</td>
1315  <td>95</td>
1316  <td>95</td>
1317  <td>88</td>
1318  <td>83</td>
1319  <td>95</td>
1320 </tr>
1321 <tr>
1322  <td>3</td>
1323  <td>37507</td>
1324  <td>107</td>
1325  <td>105</td>
1326  <td>104</td>
1327  <td>99</td>
1328  <td>100</td>
1329  <td>96</td>
1330  <td>96</td>
1331  <td>90</td>
1332 </tr>
1333 <tr>
1334  <td>4</td>
1335  <td>148739</td>
1336  <td>112</td>
1337  <td>112</td>
1338  <td>111</td>
1339  <td>111</td>
1340  <td>127</td>
1341  <td>126</td>
1342  <td>115</td>
1343  <td>117</td>
1344 </tr>
1345 <tr>
1346  <td>5</td>
1347  <td>592387</td>
1348  <td>116</td>
1349  <td>115</td>
1350  <td>114</td>
1351  <td>112</td>
1352  <td>118</td>
1353  <td>120</td>
1354  <td>131</td>
1355  <td>130</td>
1356 </tr>
1357 <tr>
1358  <td>6</td>
1359  <td>2364419</td>
1360  <td>130</td>
1361  <td>126</td>
1362  <td>120</td>
1363  <td>120</td>
1364  <td>121</td>
1365  <td>122</td>
1366  <td>121</td>
1367  <td>123</td>
1368 </tr>
1369 </table>
1370 
1371 While the PETSc results show a constant number of iterations, the iterations
1372 increase when using Trilinos. This is likely because of the different settings
1373 used for the AMG preconditioner. For performance reasons we do not allow
1374 coarsening below a couple thousand unknowns. As the coarse solver is an exact
1375 solve (we are using LU by default), a change in number of levels will
1376 influence the quality of a V-cycle. Therefore, a V-cycle is closer to an exact
1377 solver for smaller problem sizes.
1378 
1379 <a name="extensions"></a>
1380 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
1381 
1382 
1383 <a name="InvestigateTrilinositerations"></a><h4>Investigate Trilinos iterations</h4>
1384 
1385 
1386 Play with the smoothers, smoothing steps, and other properties for the
1387 Trilinos AMG to achieve an optimal preconditioner.
1388 
1389 <a name="SolvetheOseenprobleminsteadoftheStokessystem"></a><h4>Solve the Oseen problem instead of the Stokes system</h4>
1390 
1391 
1392 This change requires changing the outer solver to GMRES or BiCGStab, because
1393 the system is no longer symmetric.
1394 
1395 You can prescribe the exact flow solution as @f$b@f$ in the convective term @f$b
1396 \cdot \nabla u@f$. This should give the same solution as the original problem,
1397 if you set the right hand side to zero.
1398 
1399 <a name="Adaptiverefinement"></a><h4>Adaptive refinement</h4>
1400 
1401 
1402 So far, this tutorial program refines the mesh globally in each step.
1403 Replacing the code in StokesProblem::refine_grid() by something like
1404 @code
1405 Vector<float> estimated_error_per_cell(triangulation.n_active_cells());
1406 
1407 FEValuesExtractors::Vector velocities(0);
1408 KellyErrorEstimator<dim>::estimate(
1409  dof_handler,
1410  QGauss<dim - 1>(fe.degree + 1),
1411  std::map<types::boundary_id, const Function<dim> *>(),
1412  locally_relevant_solution,
1413  estimated_error_per_cell,
1414  fe.component_mask(velocities));
1416  triangulation, estimated_error_per_cell, 0.3, 0.0);
1417 triangulation.execute_coarsening_and_refinement();
1418 @endcode
1419 makes it simple to explore adaptive mesh refinement.
1420  *
1421  *
1422 <a name="PlainProg"></a>
1423 <h1> The plain program</h1>
1424 @include "step-55.cc"
1425 */
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
IndexSet
Definition: index_set.h:74
FEValuesExtractors
Definition: fe_values_extractors.h:82
LinearAlgebraDealII::Vector
Vector< double > Vector
Definition: generic_linear_algebra.h:43
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition: fe_update_flags.h:122
VectorTools::L2_norm
@ L2_norm
Definition: vector_tools_common.h:113
SolverCG
Definition: solver_cg.h:98
TimerOutput::Scope
Definition: timer.h:554
GridRefinement
Definition: grid_refinement.h:54
FE_Q
Definition: fe_q.h:554
dealii
Definition: namespace_dealii.h:25
DataComponentInterpretation::component_is_scalar
@ component_is_scalar
Definition: data_component_interpretation.h:55
std::sin
::VectorizedArray< Number, width > sin(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5297
Triangulation
Definition: tria.h:1109
parallel::distributed::GridRefinement::refine_and_coarsen_fixed_number
void refine_and_coarsen_fixed_number(parallel::distributed::Triangulation< dim, spacedim > &tria, const ::Vector< Number > &criteria, const double top_fraction_of_cells, const double bottom_fraction_of_cells, const types::global_cell_index max_n_cells=std::numeric_limits< types::global_cell_index >::max())
Definition: grid_refinement.cc:432
std::pow
::VectorizedArray< Number, width > pow(const ::VectorizedArray< Number, width > &, const Number p)
Definition: vectorization.h:5428
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
LAPACKSupport::V
static const char V
Definition: lapack_support.h:175
TimerOutput::summary
@ summary
Definition: timer.h:605
VectorType
IndexSet::get_view
IndexSet get_view(const size_type begin, const size_type end) const
Definition: index_set.cc:211
VectorTools::integrate_difference
void integrate_difference(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const InVector &fe_function, const Function< spacedim, typename InVector::value_type > &exact_solution, OutVector &difference, const Quadrature< dim > &q, const NormType &norm, const Function< spacedim, double > *weight=nullptr, const double exponent=2.)
MPI_Comm
VectorOperation::add
@ add
Definition: vector_operation.h:53
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
BlockDynamicSparsityPattern
Definition: block_sparsity_pattern.h:521
TimerOutput::wall_times
@ wall_times
Definition: timer.h:649
VectorTools::compute_global_error
double compute_global_error(const Triangulation< dim, spacedim > &tria, const InVector &cellwise_error, const NormType &norm, const double exponent=2.)
DataOutInterface::write_vtu_with_pvtu_record
std::string write_vtu_with_pvtu_record(const std::string &directory, const std::string &filename_without_extension, const unsigned int counter, const MPI_Comm &mpi_communicator, const unsigned int n_digits_for_counter=numbers::invalid_unsigned_int, const unsigned int n_groups=0) const
Definition: data_out_base.cc:7001
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
DoFTools::always
@ always
Definition: dof_tools.h:236
LocalIntegrators::Advection::cell_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:80
update_values
@ update_values
Shape function values.
Definition: fe_update_flags.h:78
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
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)
Table
Definition: table.h:699
FEValues< dim >
Subscriptor
Definition: subscriptor.h:62
KellyErrorEstimator
Definition: error_estimator.h:262
internal::TriangulationImplementation::n_active_cells
unsigned int n_active_cells(const internal::TriangulationImplementation::NumberCache< 1 > &c)
Definition: tria.cc:12625
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
TimerOutput
Definition: timer.h:546
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
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
LinearAlgebraDealII::BlockVector
BlockVector< double > BlockVector
Definition: generic_linear_algebra.h:48
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
FEValuesExtractors::Vector
Definition: fe_values_extractors.h:150
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5412
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
DoFTools::none
@ none
Definition: dof_tools.h:232
std::exp
::VectorizedArray< Number, width > exp(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5368
PETScWrappers::PreconditionBoomerAMG::initialize
void initialize(const MatrixBase &matrix, const AdditionalData &additional_data=AdditionalData())
Definition: petsc_precondition.cc:533
ComponentSelectFunction
Definition: function.h:560
LinearAlgebraDealII::BlockSparseMatrix
BlockSparseMatrix< double > BlockSparseMatrix
Definition: generic_linear_algebra.h:58
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
parallel::distributed::Triangulation< dim >
LinearAlgebraPETSc::MPI::PreconditionAMG
PETScWrappers::PreconditionBoomerAMG PreconditionAMG
Definition: generic_linear_algebra.h:133
VectorTools::mean
@ mean
Definition: vector_tools_common.h:79
Threads::internal::call
void call(const std::function< RT()> &function, internal::return_value< RT > &ret_val)
Definition: thread_management.h:607
types
Definition: types.h:31
DoFTools::count_dofs_per_fe_block
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandlerType &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition: dof_tools.cc:2012
QGauss
Definition: quadrature_lib.h:40
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
SmartPointer
Definition: smartpointer.h:68
BlockVector::reinit
void reinit(const unsigned int n_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
value
static const bool value
Definition: dof_tools_constraints.cc:433
AffineConstraints< double >
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition: fe_update_flags.h:129
LinearAlgebraDealII::SparseMatrix
SparseMatrix< double > SparseMatrix
Definition: generic_linear_algebra.h:53
AdaptationStrategies::Refinement::split
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
Assert
#define Assert(cond, exc)
Definition: exceptions.h:1419
Utilities::MPI::this_mpi_process
unsigned int this_mpi_process(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:128
DoFTools::extract_locally_relevant_dofs
void extract_locally_relevant_dofs(const DoFHandlerType &dof_handler, IndexSet &dof_set)
Definition: dof_tools.cc:1173
GridGenerator::hyper_cube
void hyper_cube(Triangulation< dim, spacedim > &tria, const double left=0., const double right=1., const bool colorize=false)
Utilities::MPI::n_mpi_processes
unsigned int n_mpi_processes(const MPI_Comm &mpi_communicator)
Definition: mpi.cc:117
ConditionalOStream
Definition: conditional_ostream.h:81
SparsityTools::distribute_sparsity_pattern
void distribute_sparsity_pattern(DynamicSparsityPattern &dsp, const IndexSet &locally_owned_rows, const MPI_Comm &mpi_comm, const IndexSet &locally_relevant_rows)
Definition: sparsity_tools.cc:1046
SolverMinRes
Definition: solver_minres.h:72
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
Utilities::MPI::min
T min(const T &t, const MPI_Comm &mpi_communicator)
Point< dim >
FullMatrix< double >
Utilities::MPI::all_gather
std::vector< T > all_gather(const MPI_Comm &comm, const T &object_to_send)
VectorTools::interpolate
void interpolate(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const ComponentMask &component_mask=ComponentMask())
Function
Definition: function.h:151
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
SolverControl
Definition: solver_control.h:67
VectorTools::compute_mean_value
VectorType::value_type compute_mean_value(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const Quadrature< dim > &quadrature, const VectorType &v, const unsigned int component)
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
DataOut< dim >
numbers::PI
static constexpr double PI
Definition: numbers.h:237
Vector< double >
DoFRenumbering::component_wise
void component_wise(DoFHandlerType &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition: dof_renumbering.cc:633
FESystem
Definition: fe.h:44
scalar_product
constexpr ProductType< Number, OtherNumber >::type scalar_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, OtherNumber > &t2)
Definition: symmetric_tensor.h:3749
parallel
Definition: distributed.h:416
std::cos
::VectorizedArray< Number, width > cos(const ::VectorizedArray< Number, width > &)
Definition: vectorization.h:5324
Utilities::MPI::max
T max(const T &t, const MPI_Comm &mpi_communicator)
DataComponentInterpretation::component_is_part_of_vector
@ component_is_part_of_vector
Definition: data_component_interpretation.h:61
Utilities::MPI::MPI_InitFinalize
Definition: mpi.h:828
DataOut_DoFData::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