1896 *
std::pair<unsigned int, double>
1943 * @
ref step_18
"step-18",
deal.II
's native quadrature point data manager is employed
1947 * CellDataStorage<typename Triangulation<dim>::cell_iterator,
1948 * PointHistory<dim>>
1949 * quadrature_point_history;
1953 * A description of the finite-element system including the displacement
1954 * polynomial degree, the degree-of-freedom handler, number of DoFs per
1955 * cell and the extractor objects used to retrieve information from the
1959 * const unsigned int degree;
1960 * const FESystem<dim> fe;
1961 * DoFHandler<dim> dof_handler;
1962 * const unsigned int dofs_per_cell;
1966 * Description of how the block-system is arranged. There are 3 blocks,
1967 * the first contains a vector DOF @f$\mathbf{u}@f$ while the other two
1968 * describe scalar DOFs, @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$.
1971 * static constexpr unsigned int n_blocks = 3;
1972 * static constexpr unsigned int n_components = dim + 2;
1973 * static constexpr unsigned int first_u_component = 0;
1974 * static constexpr unsigned int p_component = dim;
1975 * static constexpr unsigned int J_component = dim + 1;
1977 * static constexpr FEValuesExtractors::Vector u_fe =
1978 * FEValuesExtractors::Vector(first_u_component);
1979 * static constexpr FEValuesExtractors::Scalar p_fe =
1980 * FEValuesExtractors::Scalar(p_component);
1981 * static constexpr FEValuesExtractors::Scalar J_fe =
1982 * FEValuesExtractors::Scalar(J_component);
1991 * std::vector<types::global_dof_index> dofs_per_block;
1992 * std::vector<types::global_dof_index> element_indices_u;
1993 * std::vector<types::global_dof_index> element_indices_p;
1994 * std::vector<types::global_dof_index> element_indices_J;
1998 * Rules for Gauss-quadrature on both the cell and faces. The number of
1999 * quadrature points on both cells and faces is recorded.
2002 * const QGauss<dim> qf_cell;
2003 * const QGauss<dim - 1> qf_face;
2004 * const unsigned int n_q_points;
2005 * const unsigned int n_q_points_f;
2009 * Objects that store the converged solution and right-hand side vectors,
2010 * as well as the tangent matrix. There is an AffineConstraints object used
2011 * to keep track of constraints. We make use of a sparsity pattern
2012 * designed for a block system.
2015 * AffineConstraints<double> constraints;
2016 * BlockSparsityPattern sparsity_pattern;
2017 * BlockSparseMatrix<double> tangent_matrix;
2018 * BlockVector<double> system_rhs;
2019 * BlockVector<double> solution_n;
2023 * Then define a number of variables to store norms and update norms and
2024 * normalization factors.
2043 * void normalize(const Errors &rhs)
2045 * if (rhs.norm != 0.0)
2055 * double norm, u, p, J;
2058 * Errors error_residual, error_residual_0, error_residual_norm, error_update,
2059 * error_update_0, error_update_norm;
2063 * Methods to calculate error measures
2066 * void get_error_residual(Errors &error_residual);
2068 * void get_error_update(const BlockVector<double> &newton_update,
2069 * Errors &error_update);
2071 * std::pair<double, double> get_error_dilation() const;
2075 * Compute the volume in the spatial configuration
2078 * double compute_vol_current() const;
2082 * Print information to screen in a pleasing way...
2085 * static void print_conv_header();
2087 * void print_conv_footer();
2093 * <a name="step_44-ImplementationofthecodeSolidcodeclass"></a>
2094 * <h3>Implementation of the <code>Solid</code> class</h3>
2099 * <a name="step_44-Publicinterface"></a>
2100 * <h4>Public interface</h4>
2104 * We initialize the Solid class using data extracted from the parameter file.
2107 * template <int dim>
2108 * Solid<dim>::Solid(const std::string &input_file)
2109 * : parameters(input_file)
2110 * , vol_reference(0.)
2111 * , triangulation(Triangulation<dim>::maximum_smoothing)
2112 * , time(parameters.end_time, parameters.delta_t)
2113 * , timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
2114 * , degree(parameters.poly_degree)
2118 * The Finite Element System is composed of dim continuous displacement
2119 * DOFs, and discontinuous pressure and dilatation DOFs. In an attempt to
2120 * satisfy the Babuska-Brezzi or LBB stability conditions (see Hughes
2121 * (2000)), we set up a @f$Q_n \times DGP_{n-1} \times DGP_{n-1}@f$
2122 * system. @f$Q_2 \times DGP_1 \times DGP_1@f$ elements satisfy this
2123 * condition, while @f$Q_1 \times DGP_0 \times DGP_0@f$ elements do
2124 * not. However, it has been shown that the latter demonstrate good
2125 * convergence characteristics nonetheless.
2128 * fe(FE_Q<dim>(parameters.poly_degree) ^ dim, // displacement
2129 * FE_DGP<dim>(parameters.poly_degree - 1), // pressure
2130 * FE_DGP<dim>(parameters.poly_degree - 1)) // dilatation
2131 * , dof_handler(triangulation)
2132 * , dofs_per_cell(fe.n_dofs_per_cell())
2133 * , dofs_per_block(n_blocks)
2134 * , qf_cell(parameters.quad_order)
2135 * , qf_face(parameters.quad_order)
2136 * , n_q_points(qf_cell.size())
2137 * , n_q_points_f(qf_face.size())
2139 * Assert(dim == 2 || dim == 3,
2140 * ExcMessage("This problem only works in 2 or 3 space dimensions."));
2144 * Next we compute some information from the FE system that describes which
2145 * local element DOFs are attached to which block component. This is used
2146 * later to extract sub-blocks from the global matrix.
2150 * In essence, all we need is for the FESystem object to indicate to which
2151 * block component a DOF is attached. We can do that via the
2152 * FiniteElement::shape_function_belongs_to() function.
2155 * for (unsigned int k = 0; k < fe.n_dofs_per_cell(); ++k)
2157 * if (fe.shape_function_belongs_to(k, u_fe))
2158 * element_indices_u.push_back(k);
2159 * else if (fe.shape_function_belongs_to(k, p_fe))
2160 * element_indices_p.push_back(k);
2161 * else if (fe.shape_function_belongs_to(k, J_fe))
2162 * element_indices_J.push_back(k);
2164 * DEAL_II_ASSERT_UNREACHABLE();
2171 * In solving the quasi-static problem, the time becomes a loading parameter,
2172 * i.e. we increasing the loading linearly with time, making the two concepts
2173 * interchangeable. We choose to increment time linearly using a constant time
2178 * We start the function with preprocessing, setting the initial dilatation
2179 * values, and then output the initial grid before starting the simulation
2180 * proper with the first time (and loading)
2185 * Care must be taken (or at least some thought given) when imposing the
2186 * constraint @f$\widetilde{J}=1@f$ on the initial solution field. The constraint
2187 * corresponds to the determinant of the deformation gradient in the
2188 * undeformed configuration, which is the identity tensor. We use
2189 * FE_DGP bases to interpolate the dilatation field, thus we can't
2191 * coefficients
of a
truncated Legendre polynomial.
2205 *
template <
int dim>
2212 *
constraints.close();
2234 *
while (time.current() < time.end())
2263 * <a name=
"step_44-Privateinterface"></a>
2269 * <a name=
"step_44-Threadingbuildingblocksstructures"></a>
2294 *
std::vector<types::global_dof_index> local_dof_indices;
2299 *
, local_dof_indices(dofs_per_cell)
2318 *
template <
int dim>
2324 *
std::vector<std::vector<double>>
Nx;
2325 *
std::vector<std::vector<Tensor<2, dim>>>
grad_Nx;
2326 *
std::vector<std::vector<SymmetricTensor<2, dim>>>
symm_grad_Nx;
2344 *
: fe_values(
rhs.fe_values.get_fe(),
2345 *
rhs.fe_values.get_quadrature(),
2346 *
rhs.fe_values.get_update_flags())
2347 *
, fe_face_values(
rhs.fe_face_values.get_fe(),
2348 *
rhs.fe_face_values.get_quadrature(),
2349 *
rhs.fe_face_values.get_update_flags())
2357 *
const unsigned int n_q_points =
Nx.size();
2358 *
const unsigned int n_dofs_per_cell =
Nx[0].size();
2365 *
for (
unsigned int k = 0;
k < n_dofs_per_cell; ++
k)
2392 *
template <
int dim>
2396 *
std::vector<types::global_dof_index> local_dof_indices;
2409 *
const unsigned int n_u,
2410 *
const unsigned int n_p,
2411 *
const unsigned int n_J)
2413 *
, local_dof_indices(dofs_per_cell)
2414 *
,
k_orig(dofs_per_cell, dofs_per_cell)
2438 *
template <
int dim>
2465 *
template <
int dim>
2476 * vector
so that we don't have to copy this large data structure. We then
2477 * define a number of vectors to extract the solution values and gradients at
2478 * the quadrature points.
2481 * template <int dim>
2482 * struct Solid<dim>::ScratchData_UQPH
2484 * const BlockVector<double> &solution_total;
2486 * std::vector<Tensor<2, dim>> solution_grads_u_total;
2487 * std::vector<double> solution_values_p_total;
2488 * std::vector<double> solution_values_J_total;
2490 * FEValues<dim> fe_values;
2492 * ScratchData_UQPH(const FiniteElement<dim> &fe_cell,
2493 * const QGauss<dim> &qf_cell,
2494 * const UpdateFlags uf_cell,
2495 * const BlockVector<double> &solution_total)
2496 * : solution_total(solution_total)
2497 * , solution_grads_u_total(qf_cell.size())
2498 * , solution_values_p_total(qf_cell.size())
2499 * , solution_values_J_total(qf_cell.size())
2500 * , fe_values(fe_cell, qf_cell, uf_cell)
2503 * ScratchData_UQPH(const ScratchData_UQPH &rhs)
2504 * : solution_total(rhs.solution_total)
2505 * , solution_grads_u_total(rhs.solution_grads_u_total)
2506 * , solution_values_p_total(rhs.solution_values_p_total)
2507 * , solution_values_J_total(rhs.solution_values_J_total)
2508 * , fe_values(rhs.fe_values.get_fe(),
2509 * rhs.fe_values.get_quadrature(),
2510 * rhs.fe_values.get_update_flags())
2515 * const unsigned int n_q_points = solution_grads_u_total.size();
2516 * for (unsigned int q = 0; q < n_q_points; ++q)
2518 * solution_grads_u_total[q] = 0.0;
2519 * solution_values_p_total[q] = 0.0;
2520 * solution_values_J_total[q] = 0.0;
2529 * <a name="step_44-Solidmake_grid"></a>
2530 * <h4>Solid::make_grid</h4>
2534 * On to the first of the private member functions. Here we create the
2535 * triangulation of the domain, for which we choose the scaled cube with each
2536 * face given a boundary ID number. The grid must be refined at least once
2537 * for the indentation problem.
2541 * We then determine the volume of the reference configuration and print it
2545 * template <int dim>
2546 * void Solid<dim>::make_grid()
2548 * GridGenerator::hyper_rectangle(
2550 * (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
2551 * (dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
2553 * GridTools::scale(parameters.scale, triangulation);
2554 * triangulation.refine_global(std::max(1U, parameters.global_refinement));
2556 * vol_reference = GridTools::volume(triangulation);
2557 * std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
2561 * Since we wish to apply a Neumann BC to a patch on the top surface, we
2562 * must find the cell faces in this part of the domain and mark them with
2563 * a distinct boundary ID number. The faces we are looking for are on the
2564 * +y surface and will get boundary ID 6 (zero through five are already
2565 * used when creating the six faces of the cube domain):
2568 * for (const auto &cell : triangulation.active_cell_iterators())
2569 * for (const auto &face : cell->face_iterators())
2571 * if (face->at_boundary() == true &&
2572 * face->center()[1] == 1.0 * parameters.scale)
2576 * if (face->center()[0] < 0.5 * parameters.scale &&
2577 * face->center()[2] < 0.5 * parameters.scale)
2578 * face->set_boundary_id(6);
2582 * if (face->center()[0] < 0.5 * parameters.scale)
2583 * face->set_boundary_id(6);
2593 * <a name="step_44-Solidsystem_setup"></a>
2594 * <h4>Solid::system_setup</h4>
2598 * Next we describe how the FE system is setup. We first determine the number
2599 * of components per block. Since the displacement is a vector component, the
2600 * first dim components belong to it, while the next two describe scalar
2601 * pressure and dilatation DOFs.
2604 * template <int dim>
2605 * void Solid<dim>::system_setup()
2607 * timer.enter_subsection("Setup system");
2609 * std::vector<unsigned int> block_component(n_components,
2610 * u_dof); // Displacement
2611 * block_component[p_component] = p_dof; // Pressure
2612 * block_component[J_component] = J_dof; // Dilatation
2616 * The DOF handler is then initialized and we renumber the grid in an
2617 * efficient manner. We also record the number of DOFs per block.
2620 * dof_handler.distribute_dofs(fe);
2621 * DoFRenumbering::Cuthill_McKee(dof_handler);
2622 * DoFRenumbering::component_wise(dof_handler, block_component);
2625 * DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
2627 * std::cout << "Triangulation:"
2628 * << "\n\t Number of active cells: "
2629 * << triangulation.n_active_cells()
2630 * << "\n\t Number of degrees of freedom: " << dof_handler.n_dofs()
2635 * Setup the sparsity pattern and tangent matrix
2638 * tangent_matrix.clear();
2640 * BlockDynamicSparsityPattern dsp(dofs_per_block, dofs_per_block);
2644 * The global system matrix initially has the following structure
2646 * \underbrace{\begin{bmatrix}
2647 * \mathsf{\mathbf{K}}_{uu} & \mathsf{\mathbf{K}}_{u\widetilde{p}} &
2649 * \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
2650 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}
2651 * \\ \mathbf{0} & \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
2652 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
2653 * \end{bmatrix}}_{\mathsf{\mathbf{K}}(\mathbf{\Xi}_{\textrm{i}})}
2654 * \underbrace{\begin{bmatrix}
2656 * \\ d \widetilde{\mathsf{\mathbf{p}}}
2657 * \\ d \widetilde{\mathsf{\mathbf{J}}}
2658 * \end{bmatrix}}_{d \mathbf{\Xi}}
2660 * \underbrace{\begin{bmatrix}
2661 * \mathsf{\mathbf{F}}_{u}(\mathbf{u}_{\textrm{i}})
2662 * \\ \mathsf{\mathbf{F}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}})
2663 * \\ \mathsf{\mathbf{F}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}})
2664 * \end{bmatrix}}_{ \mathsf{\mathbf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, .
2666 * We optimize the sparsity pattern to reflect this structure
2667 * and prevent unnecessary data creation for the right-diagonal
2671 * Table<2, DoFTools::Coupling> coupling(n_components, n_components);
2672 * for (unsigned int ii = 0; ii < n_components; ++ii)
2673 * for (unsigned int jj = 0; jj < n_components; ++jj)
2674 * if (((ii < p_component) && (jj == J_component)) ||
2675 * ((ii == J_component) && (jj < p_component)) ||
2676 * ((ii == p_component) && (jj == p_component)))
2677 * coupling[ii][jj] = DoFTools::none;
2679 * coupling[ii][jj] = DoFTools::always;
2680 * DoFTools::make_sparsity_pattern(
2681 * dof_handler, coupling, dsp, constraints, false);
2682 * sparsity_pattern.copy_from(dsp);
2685 * tangent_matrix.reinit(sparsity_pattern);
2689 * We then set up storage vectors
2692 * system_rhs.reinit(dofs_per_block);
2693 * solution_n.reinit(dofs_per_block);
2697 * ...and finally set up the quadrature
2703 * timer.leave_subsection();
2710 * <a name="step_44-Solidsetup_qph"></a>
2711 * <h4>Solid::setup_qph</h4>
2712 * The method used to store quadrature information is already described in
2713 * @ref step_18 "step-18". Here we implement a similar setup for a SMP machine.
2717 * Firstly the actual QPH data objects are created. This must be done only
2718 * once the grid is refined to its finest level.
2721 * template <int dim>
2722 * void Solid<dim>::setup_qph()
2724 * std::cout << " Setting up quadrature point data..." << std::endl;
2726 * quadrature_point_history.initialize(triangulation.begin_active(),
2727 * triangulation.end(),
2732 * Next we set up the initial quadrature point data.
2733 * Note that when the quadrature point data is retrieved,
2734 * it is returned as a vector of smart pointers.
2737 * for (const auto &cell : triangulation.active_cell_iterators())
2739 * const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
2740 * quadrature_point_history.get_data(cell);
2741 * AssertDimension(lqph.size(), n_q_points);
2743 * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2744 * lqph[q_point]->setup_lqp(parameters);
2751 * <a name="step_44-Solidupdate_qph_incremental"></a>
2752 * <h4>Solid::update_qph_incremental</h4>
2753 * As the update of QP information occurs frequently and involves a number of
2754 * expensive operations, we define a multithreaded approach to distributing
2755 * the task across a number of CPU cores.
2759 * To start this, we first we need to obtain the total solution as it stands
2760 * at this Newton increment and then create the initial copy of the scratch
2761 * and copy data objects:
2764 * template <int dim>
2766 * Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
2768 * timer.enter_subsection("Update QPH data");
2769 * std::cout << " UQPH " << std::flush;
2771 * const BlockVector<double> solution_total(
2772 * get_total_solution(solution_delta));
2774 * const UpdateFlags uf_UQPH(update_values | update_gradients);
2775 * PerTaskData_UQPH per_task_data_UQPH;
2776 * ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
2780 * We then pass them and the one-cell update function to the WorkStream to
2784 * WorkStream::run(dof_handler.active_cell_iterators(),
2786 * &Solid::update_qph_incremental_one_cell,
2787 * &Solid::copy_local_to_global_UQPH,
2788 * scratch_data_UQPH,
2789 * per_task_data_UQPH);
2791 * timer.leave_subsection();
2797 * Now we describe how we extract data from the solution vector and pass it
2798 * along to each QP storage object for processing.
2801 * template <int dim>
2802 * void Solid<dim>::update_qph_incremental_one_cell(
2803 * const typename DoFHandler<dim>::active_cell_iterator &cell,
2804 * ScratchData_UQPH &scratch,
2805 * PerTaskData_UQPH & /*data*/)
2807 * const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
2808 * quadrature_point_history.get_data(cell);
2809 * AssertDimension(lqph.size(), n_q_points);
2811 * AssertDimension(scratch.solution_grads_u_total.size(), n_q_points);
2812 * AssertDimension(scratch.solution_values_p_total.size(), n_q_points);
2813 * AssertDimension(scratch.solution_values_J_total.size(), n_q_points);
2819 * We first need to find the values and gradients at quadrature points
2820 * inside the current cell and then we update each local QP using the
2821 * displacement gradient and total pressure and dilatation solution
2825 * scratch.fe_values.reinit(cell);
2826 * scratch.fe_values[u_fe].get_function_gradients(
2827 * scratch.solution_total, scratch.solution_grads_u_total);
2828 * scratch.fe_values[p_fe].get_function_values(
2829 * scratch.solution_total, scratch.solution_values_p_total);
2830 * scratch.fe_values[J_fe].get_function_values(
2831 * scratch.solution_total, scratch.solution_values_J_total);
2833 * for (const unsigned int q_point :
2834 * scratch.fe_values.quadrature_point_indices())
2835 * lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point],
2836 * scratch.solution_values_p_total[q_point],
2837 * scratch.solution_values_J_total[q_point]);
2844 * <a name="step_44-Solidsolve_nonlinear_timestep"></a>
2845 * <h4>Solid::solve_nonlinear_timestep</h4>
2849 * The next function is the driver method for the Newton-Raphson scheme. At
2850 * its top we create a new vector to store the current Newton update step,
2851 * reset the error storage objects and print solver header.
2854 * template <int dim>
2855 * void Solid<dim>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
2857 * std::cout << std::endl
2858 * << "Timestep " << time.get_timestep() << " @ " << time.current()
2859 * << 's
' << std::endl;
2861 * BlockVector<double> newton_update(dofs_per_block);
2863 * error_residual.reset();
2864 * error_residual_0.reset();
2865 * error_residual_norm.reset();
2866 * error_update.reset();
2867 * error_update_0.reset();
2868 * error_update_norm.reset();
2870 * print_conv_header();
2874 * We now perform a number of Newton iterations to iteratively solve the
2875 * nonlinear problem. Since the problem is fully nonlinear and we are
2876 * using a full Newton method, the data stored in the tangent matrix and
2877 * right-hand side vector is not reusable and must be cleared at each
2878 * Newton step. We then initially build the linear system and
2879 * check for convergence (and store this value in the first iteration).
2880 * The unconstrained DOFs of the rhs vector hold the out-of-balance
2881 * forces, and collectively determine whether or not the equilibrium
2882 * solution has been attained.
2886 * Although for this particular problem we could potentially construct the
2887 * RHS vector before assembling the system matrix, for the sake of
2888 * extensibility we choose not to do so. The benefit to assembling the RHS
2889 * vector and system matrix separately is that the latter is an expensive
2890 * operation and we can potentially avoid an extra assembly process by not
2891 * assembling the tangent matrix when convergence is attained. However, this
2892 * makes parallelizing the code using MPI more difficult. Furthermore, when
2893 * extending the problem to the transient case additional contributions to
2894 * the RHS may result from the time discretization and application of
2895 * constraints for the velocity and acceleration fields.
2898 * unsigned int newton_iteration = 0;
2899 * for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration)
2901 * std::cout << ' ' << std::setw(2) << newton_iteration << ' '
2906 * We construct the linear system, but hold off on solving it
2907 * (a step that should be significantly more expensive than assembly):
2910 * make_constraints(newton_iteration);
2911 * assemble_system();
2915 * We can now determine the normalized residual error and check for
2916 * solution convergence:
2919 * get_error_residual(error_residual);
2920 * if (newton_iteration == 0)
2921 * error_residual_0 = error_residual;
2923 * error_residual_norm = error_residual;
2924 * error_residual_norm.normalize(error_residual_0);
2926 * if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u &&
2927 * error_residual_norm.u <= parameters.tol_f)
2929 * std::cout << " CONVERGED! " << std::endl;
2930 * print_conv_footer();
2937 * If we have decided that we want to continue with the iteration, we
2938 * solve the linearized system:
2941 * const std::pair<unsigned int, double> lin_solver_output =
2942 * solve_linear_system(newton_update);
2946 * We can now determine the normalized Newton update error:
2949 * get_error_update(newton_update, error_update);
2950 * if (newton_iteration == 0)
2951 * error_update_0 = error_update;
2953 * error_update_norm = error_update;
2954 * error_update_norm.normalize(error_update_0);
2958 * Lastly, since we implicitly accept the solution step we can perform
2959 * the actual update of the solution increment for the current time
2960 * step, update all quadrature point information pertaining to
2961 * this new displacement and stress state and continue iterating:
2964 * solution_delta += newton_update;
2965 * update_qph_incremental(solution_delta);
2967 * std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
2968 * << std::scientific << lin_solver_output.first << " "
2969 * << lin_solver_output.second << " "
2970 * << error_residual_norm.norm << " " << error_residual_norm.u
2971 * << " " << error_residual_norm.p << " "
2972 * << error_residual_norm.J << " " << error_update_norm.norm
2973 * << " " << error_update_norm.u << " " << error_update_norm.p
2974 * << " " << error_update_norm.J << " " << std::endl;
2979 * At the end, if it turns out that we have in fact done more iterations
2980 * than the parameter file allowed, we raise an exception that can be
2981 * caught in the main() function. The call <code>AssertThrow(condition,
2982 * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
2983 * exc_object;</code> but the former form fills certain fields in the
2984 * exception object that identify the location (filename and line number)
2985 * where the exception was raised to make it simpler to identify where the
2989 * AssertThrow(newton_iteration < parameters.max_iterations_NR,
2990 * ExcMessage("No convergence in nonlinear solver!"));
2997 * <a name="step_44-Solidprint_conv_headerandSolidprint_conv_footer"></a>
2998 * <h4>Solid::print_conv_header and Solid::print_conv_footer</h4>
3002 * This program prints out data in a nice table that is updated
3003 * on a per-iteration basis. The next two functions set up the table
3004 * header and footer:
3007 * template <int dim>
3008 * void Solid<dim>::print_conv_header()
3010 * static const unsigned int l_width = 150;
3012 * for (unsigned int i = 0; i < l_width; ++i)
3014 * std::cout << std::endl;
3016 * std::cout << " SOLVER STEP "
3017 * << " | LIN_IT LIN_RES RES_NORM "
3018 * << " RES_U RES_P RES_J NU_NORM "
3019 * << " NU_U NU_P NU_J " << std::endl;
3021 * for (unsigned int i = 0; i < l_width; ++i)
3023 * std::cout << std::endl;
3028 * template <int dim>
3029 * void Solid<dim>::print_conv_footer()
3031 * static const unsigned int l_width = 150;
3033 * for (unsigned int i = 0; i < l_width; ++i)
3035 * std::cout << std::endl;
3037 * const std::pair<double, double> error_dil = get_error_dilation();
3039 * std::cout << "Relative errors:" << std::endl
3040 * << "Displacement:\t" << error_update.u / error_update_0.u
3042 * << "Force: \t\t" << error_residual.u / error_residual_0.u
3044 * << "Dilatation:\t" << error_dil.first << std::endl
3045 * << "v / V_0:\t" << error_dil.second * vol_reference << " / "
3046 * << vol_reference << " = " << error_dil.second << std::endl;
3053 * <a name="step_44-Solidget_error_dilation"></a>
3054 * <h4>Solid::get_error_dilation</h4>
3058 * Calculate the volume of the domain in the spatial configuration
3061 * template <int dim>
3062 * double Solid<dim>::compute_vol_current() const
3064 * double vol_current = 0.0;
3066 * FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
3068 * for (const auto &cell : triangulation.active_cell_iterators())
3070 * fe_values.reinit(cell);
3074 * In contrast to that which was previously called for,
3075 * in this instance the quadrature point data is specifically
3076 * non-modifiable since we will only be accessing data.
3077 * We ensure that the right get_data function is called by
3078 * marking this update function as constant.
3081 * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3082 * quadrature_point_history.get_data(cell);
3083 * AssertDimension(lqph.size(), n_q_points);
3085 * for (const unsigned int q_point : fe_values.quadrature_point_indices())
3087 * const double det_F_qp = lqph[q_point]->get_det_F();
3088 * const double JxW = fe_values.JxW(q_point);
3090 * vol_current += det_F_qp * JxW;
3093 * Assert(vol_current > 0.0, ExcInternalError());
3094 * return vol_current;
3099 * Calculate how well the dilatation @f$\widetilde{J}@f$ agrees with @f$J
3100 * \dealcoloneq \textrm{det}\ \mathbf{F}@f$ from the @f$L^2@f$ error @f$ \bigl[
3101 * \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}@f$.
3102 * We also return the ratio of the current volume of the
3103 * domain to the reference volume. This is of interest for incompressible
3104 * media where we want to check how well the isochoric constraint has been
3108 * template <int dim>
3109 * std::pair<double, double> Solid<dim>::get_error_dilation() const
3111 * double dil_L2_error = 0.0;
3113 * FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
3115 * for (const auto &cell : triangulation.active_cell_iterators())
3117 * fe_values.reinit(cell);
3119 * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3120 * quadrature_point_history.get_data(cell);
3121 * AssertDimension(lqph.size(), n_q_points);
3123 * for (const unsigned int q_point : fe_values.quadrature_point_indices())
3125 * const double det_F_qp = lqph[q_point]->get_det_F();
3126 * const double J_tilde_qp = lqph[q_point]->get_J_tilde();
3127 * const double the_error_qp_squared =
3128 * Utilities::fixed_power<2>((det_F_qp - J_tilde_qp));
3129 * const double JxW = fe_values.JxW(q_point);
3131 * dil_L2_error += the_error_qp_squared * JxW;
3135 * return std::make_pair(std::sqrt(dil_L2_error),
3136 * compute_vol_current() / vol_reference);
3143 * <a name="step_44-Solidget_error_residual"></a>
3144 * <h4>Solid::get_error_residual</h4>
3148 * Determine the true residual error for the problem. That is, determine the
3149 * error in the residual for the unconstrained degrees of freedom. Note that
3150 * to do so, we need to ignore constrained DOFs by setting the residual in
3151 * these vector components to zero.
3154 * template <int dim>
3155 * void Solid<dim>::get_error_residual(Errors &error_residual)
3157 * BlockVector<double> error_res(dofs_per_block);
3159 * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
3160 * if (!constraints.is_constrained(i))
3161 * error_res(i) = system_rhs(i);
3163 * error_residual.norm = error_res.l2_norm();
3164 * error_residual.u = error_res.block(u_dof).l2_norm();
3165 * error_residual.p = error_res.block(p_dof).l2_norm();
3166 * error_residual.J = error_res.block(J_dof).l2_norm();
3173 * <a name="step_44-Solidget_error_update"></a>
3174 * <h4>Solid::get_error_update</h4>
3178 * Determine the true Newton update error for the problem
3181 * template <int dim>
3182 * void Solid<dim>::get_error_update(const BlockVector<double> &newton_update,
3183 * Errors &error_update)
3185 * BlockVector<double> error_ud(dofs_per_block);
3186 * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
3187 * if (!constraints.is_constrained(i))
3188 * error_ud(i) = newton_update(i);
3190 * error_update.norm = error_ud.l2_norm();
3191 * error_update.u = error_ud.block(u_dof).l2_norm();
3192 * error_update.p = error_ud.block(p_dof).l2_norm();
3193 * error_update.J = error_ud.block(J_dof).l2_norm();
3201 * <a name="step_44-Solidget_total_solution"></a>
3202 * <h4>Solid::get_total_solution</h4>
3206 * This function provides the total solution, which is valid at any Newton
3207 * step. This is required as, to reduce computational error, the total
3208 * solution is only updated at the end of the timestep.
3211 * template <int dim>
3212 * BlockVector<double> Solid<dim>::get_total_solution(
3213 * const BlockVector<double> &solution_delta) const
3215 * BlockVector<double> solution_total(solution_n);
3216 * solution_total += solution_delta;
3217 * return solution_total;
3224 * <a name="step_44-Solidassemble_system"></a>
3225 * <h4>Solid::assemble_system</h4>
3229 * Since we use TBB for assembly, we simply setup a copy of the
3230 * data structures required for the process and pass them, along
3231 * with the assembly functions to the WorkStream object for processing. Note
3232 * that we must ensure that the matrix and RHS vector are reset before any
3233 * assembly operations can occur. Furthermore, since we are describing a
3234 * problem with Neumann BCs, we will need the face normals and so must specify
3235 * this in the face update flags.
3238 * template <int dim>
3239 * void Solid<dim>::assemble_system()
3241 * timer.enter_subsection("Assemble system");
3242 * std::cout << " ASM_SYS " << std::flush;
3244 * tangent_matrix = 0.0;
3247 * const UpdateFlags uf_cell(update_values | update_gradients |
3248 * update_JxW_values);
3249 * const UpdateFlags uf_face(update_values | update_normal_vectors |
3250 * update_JxW_values);
3252 * PerTaskData_ASM per_task_data(dofs_per_cell);
3253 * ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
3257 * The syntax used here to pass data to the WorkStream class
3258 * is discussed in @ref step_13 "step-13".
3262 * dof_handler.active_cell_iterators(),
3263 * [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
3264 * ScratchData_ASM &scratch,
3265 * PerTaskData_ASM &data) {
3266 * this->assemble_system_one_cell(cell, scratch, data);
3268 * [this](const PerTaskData_ASM &data) {
3269 * this->constraints.distribute_local_to_global(data.cell_matrix,
3271 * data.local_dof_indices,
3278 * timer.leave_subsection();
3283 * Of course, we still have to define how we assemble the tangent matrix
3284 * contribution for a single cell. We first need to reset and initialize some
3285 * of the scratch data structures and retrieve some basic information
3286 * regarding the DOF numbering on this cell. We can precalculate the cell
3287 * shape function values and gradients. Note that the shape function gradients
3288 * are defined with regard to the current configuration. That is
3289 * @f$\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi}
3290 * \ \mathbf{F}^{-1}@f$.
3293 * template <int dim>
3294 * void Solid<dim>::assemble_system_one_cell(
3295 * const typename DoFHandler<dim>::active_cell_iterator &cell,
3296 * ScratchData_ASM &scratch,
3297 * PerTaskData_ASM &data) const
3301 * scratch.fe_values.reinit(cell);
3302 * cell->get_dof_indices(data.local_dof_indices);
3304 * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3305 * quadrature_point_history.get_data(cell);
3306 * AssertDimension(lqph.size(), n_q_points);
3308 * for (const unsigned int q_point :
3309 * scratch.fe_values.quadrature_point_indices())
3311 * const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
3312 * for (const unsigned int k : scratch.fe_values.dof_indices())
3314 * const unsigned int k_group = fe.system_to_base_index(k).first.first;
3316 * if (k_group == u_dof)
3318 * scratch.grad_Nx[q_point][k] =
3319 * scratch.fe_values[u_fe].gradient(k, q_point) * F_inv;
3320 * scratch.symm_grad_Nx[q_point][k] =
3321 * symmetrize(scratch.grad_Nx[q_point][k]);
3323 * else if (k_group == p_dof)
3324 * scratch.Nx[q_point][k] =
3325 * scratch.fe_values[p_fe].value(k, q_point);
3326 * else if (k_group == J_dof)
3327 * scratch.Nx[q_point][k] =
3328 * scratch.fe_values[J_fe].value(k, q_point);
3330 * DEAL_II_ASSERT_UNREACHABLE();
3336 * Now we build the local cell @ref GlossStiffnessMatrix "stiffness matrix" and RHS vector. Since the
3337 * global and local system matrices are symmetric, we can exploit this
3338 * property by building only the lower half of the local matrix and copying
3339 * the values to the upper half. So we only assemble half of the
3340 * @f$\mathsf{\mathbf{k}}_{uu}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
3341 * \widetilde{p}} = \mathbf{0}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{J}
3342 * \widetilde{J}}@f$ blocks, while the whole
3343 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$,
3344 * @f$\mathsf{\mathbf{k}}_{u \widetilde{J}} = \mathbf{0}@f$,
3345 * @f$\mathsf{\mathbf{k}}_{u \widetilde{p}}@f$ blocks are built.
3349 * In doing so, we first extract some configuration dependent variables
3350 * from our quadrature history objects for the current quadrature point.
3353 * for (const unsigned int q_point :
3354 * scratch.fe_values.quadrature_point_indices())
3356 * const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau();
3357 * const Tensor<2, dim> tau_ns = lqph[q_point]->get_tau();
3358 * const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc();
3359 * const double det_F = lqph[q_point]->get_det_F();
3360 * const double p_tilde = lqph[q_point]->get_p_tilde();
3361 * const double J_tilde = lqph[q_point]->get_J_tilde();
3362 * const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ();
3363 * const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2();
3364 * const SymmetricTensor<2, dim> &I =
3365 * Physics::Elasticity::StandardTensors<dim>::I;
3369 * These two tensors store some precomputed data. Their use will
3370 * explained shortly.
3373 * SymmetricTensor<2, dim> symm_grad_Nx_i_x_Jc;
3374 * Tensor<1, dim> grad_Nx_i_comp_i_x_tau;
3378 * Next we define some aliases to make the assembly process easier to
3382 * const std::vector<double> &N = scratch.Nx[q_point];
3383 * const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
3384 * scratch.symm_grad_Nx[q_point];
3385 * const std::vector<Tensor<2, dim>> &grad_Nx = scratch.grad_Nx[q_point];
3386 * const double JxW = scratch.fe_values.JxW(q_point);
3388 * for (const unsigned int i : scratch.fe_values.dof_indices())
3390 * const unsigned int component_i =
3391 * fe.system_to_component_index(i).first;
3392 * const unsigned int i_group = fe.system_to_base_index(i).first.first;
3396 * We first compute the contributions
3397 * from the internal forces. Note, by
3398 * definition of the rhs as the negative
3399 * of the residual, these contributions
3403 * if (i_group == u_dof)
3404 * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
3405 * else if (i_group == p_dof)
3406 * data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
3407 * else if (i_group == J_dof)
3408 * data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW;
3410 * DEAL_II_ASSERT_UNREACHABLE();
3414 * Before we go into the inner loop, we have one final chance to
3418 *
We won't be excessive here, but will rather focus on expensive
3419 * operations, namely those involving the rank-4 material stiffness
3420 * tensor and the rank-2 stress tensor.
3424 * What we may observe is that both of these tensors are contracted
3425 * with shape function gradients indexed on the "i" DoF. This
3426 * implies that this particular operation remains constant as we
3427 * loop over the "j" DoF. For that reason, we can extract this from
3428 * the inner loop and save the many operations that, for each
3429 * quadrature point and DoF index "i" and repeated over index "j"
3430 * are required to double contract a rank-2 symmetric tensor with a
3431 * rank-4 symmetric tensor, and a rank-1 tensor with a rank-2
3436 * At the loss of some readability, this small change will reduce
3437 * the assembly time of the symmetrized system by about half when
3438 * using the simulation default parameters, and becomes more
3439 * significant as the h-refinement level increases.
3442 * if (i_group == u_dof)
3444 * symm_grad_Nx_i_x_Jc = symm_grad_Nx[i] * Jc;
3445 * grad_Nx_i_comp_i_x_tau = grad_Nx[i][component_i] * tau_ns;
3453 *
for (
const unsigned int j :
3454 *
scratch.fe_values.dof_indices_ending_at(i))
3457 *
fe.system_to_component_index(
j).
first;
3458 *
const unsigned int j_group =
3459 *
fe.system_to_base_index(
j).first.first;
3485 *
data.cell_matrix(i,
j) +=
3507 *
data.cell_matrix(i,
j) -=
N[i] *
N[
j] * JxW;
3527 *
for (
const auto &face : cell->face_iterators())
3530 *
scratch.fe_face_values.
reinit(cell, face);
3533 *
scratch.fe_face_values.quadrature_point_indices())
3536 *
scratch.fe_face_values.normal_vector(
f_q_point);
3557 *
static const double p0 =
3558 *
-4.0 / (parameters.scale * parameters.scale);
3559 *
const double time_ramp = (time.current() / time.end());
3563 *
for (
const unsigned int i : scratch.fe_values.dof_indices())
3571 *
fe.system_to_component_index(i).first;
3573 *
scratch.fe_face_values.shape_value(i,
f_q_point);
3574 *
const double JxW = scratch.fe_face_values.JxW(
f_q_point);
3588 *
for (
const unsigned int i : scratch.fe_values.dof_indices())
3590 *
scratch.fe_values.dof_indices_starting_at(i + 1))
3599 * <a name=
"step_44-Solidmake_constraints"></a>
3604 * build non-homogeneous constraints in the zeroth iteration (that is, when
3605 * `apply_dirichlet_bc == true` in the code block that follows) and build
3606 * only the corresponding homogeneous constraints in the following step. While
3607 * the current example has only homogeneous constraints, previous experiences
3608 * have shown that a common error is forgetting to add the extra condition
3609 * when refactoring the code to specific uses. This could lead to errors that
3610 * are hard to debug. In this spirit, we choose to make the code more verbose
3611 * in terms of what operations are performed at each Newton step.
3614 * template <int dim>
3615 * void Solid<dim>::make_constraints(const unsigned int it_nr)
3619 * Since we (a) are dealing with an iterative Newton method, (b) are using
3620 * an incremental formulation for the displacement, and (c) apply the
3621 * constraints to the incremental displacement field, any non-homogeneous
3622 * constraints on the displacement update should only be specified at the
3623 * zeroth iteration. No subsequent contributions are to be made since the
3624 * constraints will be exactly satisfied after that iteration.
3627 * const bool apply_dirichlet_bc = (it_nr == 0);
3631 * Furthermore, after the first Newton iteration within a timestep, the
3632 * constraints remain the same and we do not need to modify or rebuild them
3633 * so long as we do not clear the @p constraints object.
3638 * std::cout << " --- " << std::flush;
3642 * std::cout << " CST " << std::flush;
3644 * if (apply_dirichlet_bc)
3648 * At the zeroth Newton iteration we wish to apply the full set of
3649 * non-homogeneous and homogeneous constraints that represent the
3650 * boundary conditions on the displacement increment. Since in general
3651 * the constraints may be different at each time step, we need to clear
3652 * the constraints matrix and completely rebuild it. An example case
3653 * would be if a surface is accelerating; in such a scenario the change
3654 * in displacement is non-constant between each time step.
3657 * constraints.clear();
3661 * The boundary conditions for the indentation problem in 3d are as
3662 * follows: On the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry
3663 * condition to allow only planar movement while the +x and +z faces
3664 * (IDs 1,5) are traction free. In this contrived problem, part of the
3665 * +y face (ID 3) is set to have no motion in the x- and z-component.
3666 * Finally, as described earlier, the other part of the +y face has an
3667 * the applied pressure but is also constrained in the x- and
3672 * In the following, we will have to tell the function interpolation
3673 * boundary values which components of the solution vector should be
3790 *
if (constraints.has_inhomogeneities())
3798 * back
to the main @p constraints
object.
3802 *
for (
unsigned int dof = 0; dof != dof_handler.n_dofs(); ++dof)
3806 *
constraints.
clear();
3811 *
constraints.close();
3817 * <a name=
"step_44-Solidassemble_sc"></a>
3847 *
from each element
's contributions. These contributions are then added to
3848 * the global stiffness matrix. Given this description, the following two
3849 * functions should be clear:
3852 * template <int dim>
3853 * void Solid<dim>::assemble_sc()
3855 * timer.enter_subsection("Perform static condensation");
3856 * std::cout << " ASM_SC " << std::flush;
3858 * PerTaskData_SC per_task_data(dofs_per_cell,
3859 * element_indices_u.size(),
3860 * element_indices_p.size(),
3861 * element_indices_J.size());
3862 * ScratchData_SC scratch_data;
3864 * WorkStream::run(dof_handler.active_cell_iterators(),
3866 * &Solid::assemble_sc_one_cell,
3867 * &Solid::copy_local_to_global_sc,
3871 * timer.leave_subsection();
3875 * template <int dim>
3876 * void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC &data)
3878 * for (unsigned int i = 0; i < dofs_per_cell; ++i)
3879 * for (unsigned int j = 0; j < dofs_per_cell; ++j)
3880 * tangent_matrix.add(data.local_dof_indices[i],
3881 * data.local_dof_indices[j],
3882 * data.cell_matrix(i, j));
3888 * Now we describe the static condensation process. As per usual, we must
3889 * first find out which global numbers the degrees of freedom on this cell
3890 * have and reset some data structures:
3893 * template <int dim>
3894 * void Solid<dim>::assemble_sc_one_cell(
3895 * const typename DoFHandler<dim>::active_cell_iterator &cell,
3896 * ScratchData_SC &scratch,
3897 * PerTaskData_SC &data)
3901 * cell->get_dof_indices(data.local_dof_indices);
3905 * We now extract the contribution of the dofs associated with the current
3906 * cell to the global stiffness matrix. The discontinuous nature of the
3907 * @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$ interpolations mean that their is
3908 * no coupling of the local contributions at the global level. This is not
3909 * the case with the @f$\mathbf{u}@f$ dof. In other words,
3910 * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$,
3911 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{p}}@f$ and
3912 * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$, when extracted
3913 * from the global stiffness matrix are the element contributions. This
3914 * is not the case for @f$\mathsf{\mathbf{k}}_{uu}@f$.
3918 * Note: A lower-case symbol is used to denote element stiffness matrices.
3922 * Currently the matrix corresponding to
3923 * the dof associated with the current element
3924 * (denoted somewhat loosely as @f$\mathsf{\mathbf{k}}@f$)
3928 * \mathsf{\mathbf{k}}_{uu} & \mathsf{\mathbf{k}}_{u\widetilde{p}}
3930 * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
3931 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}
3932 * \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
3933 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
3938 * We now need to modify it such that it appear as
3941 * \mathsf{\mathbf{k}}_{\textrm{con}} &
3942 * \mathsf{\mathbf{k}}_{u\widetilde{p}} & \mathbf{0}
3943 * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
3944 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
3945 * \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
3946 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
3948 * with @f$\mathsf{\mathbf{k}}_{\textrm{con}} = \bigl[
3949 * \mathsf{\mathbf{k}}_{uu} +\overline{\overline{\mathsf{\mathbf{k}}}}~
3950 * \bigr]@f$ where @f$ \overline{\overline{\mathsf{\mathbf{k}}}}
3951 * \dealcoloneq \mathsf{\mathbf{k}}_{u\widetilde{p}}
3952 * \overline{\mathsf{\mathbf{k}}} \mathsf{\mathbf{k}}_{\widetilde{p}u}
3956 * \overline{\mathsf{\mathbf{k}}} =
3957 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}}^{-1}
3958 * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}}
3959 * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
3964 * At this point, we need to take note of
3965 * the fact that global data already exists
3966 * in the @f$\mathsf{\mathbf{K}}_{uu}@f$,
3967 * @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
3969 * @f$\mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{p}}@f$
3970 * sub-blocks. So if we are to modify them, we must account for the data
3971 * that is already there (i.e. simply add to it or remove it if
3972 * necessary). Since the copy_local_to_global operation is a "+="
3973 * operation, we need to take this into account
3977 * For the @f$\mathsf{\mathbf{K}}_{uu}@f$ block in particular, this means that
3978 * contributions have been added from the surrounding cells, so we need to
3979 * be careful when we manipulate this block. We can't
just erase
the
3989 *
Since we don't have access to @f$\mathsf{\mathbf{k}}_{uu}@f$,
3990 * but we know its contribution is added to
3991 * the global @f$\mathsf{\mathbf{K}}_{uu}@f$ matrix, we just want
3992 * to add the element wise
3993 * static-condensation @f$\overline{\overline{\mathsf{\mathbf{k}}}}@f$.
3997 * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$:
3998 * Similarly, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
3999 * \widetilde{J}}@f$ exists in
4000 * the subblock. Since the copy
4001 * operation is a += operation, we
4002 * need to subtract the existing
4003 * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
4004 * submatrix in addition to
4005 * "adding" that which we wish to
4010 * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}@f$:
4011 * Since the global matrix
4012 * is symmetric, this block is the
4013 * same as the one above and we
4015 * @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
4016 * as a substitute for this one.
4020 * We first extract element data from the
4021 * system matrix. So first we get the
4022 * entire subblock for the cell, then
4023 * extract @f$\mathsf{\mathbf{k}}@f$
4024 * for the dofs associated with
4025 * the current element
4028 * data.k_orig.extract_submatrix_from(tangent_matrix,
4029 * data.local_dof_indices,
4030 * data.local_dof_indices);
4033 * and next the local matrices for
4034 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} u}@f$
4035 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
4037 * @f$\mathsf{\mathbf{k}}_{ \widetilde{J} \widetilde{J}}@f$:
4040 * data.k_pu.extract_submatrix_from(data.k_orig,
4041 * element_indices_p,
4042 * element_indices_u);
4043 * data.k_pJ.extract_submatrix_from(data.k_orig,
4044 * element_indices_p,
4045 * element_indices_J);
4046 * data.k_JJ.extract_submatrix_from(data.k_orig,
4047 * element_indices_J,
4048 * element_indices_J);
4052 * To get the inverse of @f$\mathsf{\mathbf{k}}_{\widetilde{p}
4053 * \widetilde{J}}@f$, we invert it directly. This operation is relatively
4054 * inexpensive since @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
4055 * since block-diagonal.
4058 * data.k_pJ_inv.invert(data.k_pJ);
4062 * Now we can make condensation terms to
4063 * add to the @f$\mathsf{\mathbf{k}}_{uu}@f$
4064 * block and put them in
4065 * the cell local matrix
4067 * \mathsf{\mathbf{A}}
4069 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4070 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4074 * data.k_pJ_inv.mmult(data.A, data.k_pu);
4078 * \mathsf{\mathbf{B}}
4080 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4081 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4082 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4086 * data.k_JJ.mmult(data.B, data.A);
4090 * \mathsf{\mathbf{C}}
4092 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
4093 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4094 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4095 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4099 * data.k_pJ_inv.Tmmult(data.C, data.B);
4103 * \overline{\overline{\mathsf{\mathbf{k}}}}
4105 * \mathsf{\mathbf{k}}_{u \widetilde{p}}
4106 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
4107 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4108 * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4109 * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4113 * data.k_pu.Tmmult(data.k_bbar, data.C);
4114 * data.k_bbar.scatter_matrix_to(element_indices_u,
4115 * element_indices_u,
4116 * data.cell_matrix);
4121 * @f$\mathsf{\mathbf{k}}^{-1}_{ \widetilde{p} \widetilde{J}}@f$
4123 * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
4124 * block for post-processing. Note again
4125 * that we need to remove the
4126 * contribution that already exists there.
4129 * data.k_pJ_inv.add(-1.0, data.k_pJ);
4130 * data.k_pJ_inv.scatter_matrix_to(element_indices_p,
4131 * element_indices_J,
4132 * data.cell_matrix);
4138 * <a name="step_44-Solidsolve_linear_system"></a>
4139 * <h4>Solid::solve_linear_system</h4>
4140 * We now have all of the necessary components to use one of two possible
4141 * methods to solve the linearised system. The first is to perform static
4142 * condensation on an element level, which requires some alterations
4143 * to the tangent matrix and RHS vector. Alternatively, the full block
4144 * system can be solved by performing condensation on a global level.
4145 * Below we implement both approaches.
4148 * template <int dim>
4149 * std::pair<unsigned int, double>
4150 * Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
4152 * unsigned int lin_it = 0;
4153 * double lin_res = 0.0;
4155 * if (parameters.use_static_condensation == true)
4159 * Firstly, here is the approach using the (permanent) augmentation of
4160 * the tangent matrix. For the following, recall that
4162 * \mathsf{\mathbf{K}}_{\textrm{store}}
4165 * \mathsf{\mathbf{K}}_{\textrm{con}} &
4166 * \mathsf{\mathbf{K}}_{u\widetilde{p}} & \mathbf{0}
4167 * \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
4168 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
4170 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
4171 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, .
4175 * d \widetilde{\mathsf{\mathbf{p}}}
4177 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4179 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4181 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4182 * d \widetilde{\mathsf{\mathbf{J}}} \bigr]
4183 * \\ d \widetilde{\mathsf{\mathbf{J}}}
4185 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
4187 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4188 * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
4189 * \mathsf{\mathbf{u}} \bigr]
4190 * \\ \Rightarrow d \widetilde{\mathsf{\mathbf{p}}}
4191 * &= \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4192 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4194 * \underbrace{\bigl[\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4195 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4196 * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}\bigr]}_{\overline{\mathsf{\mathbf{K}}}}\bigl[
4197 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4198 * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
4199 * \mathsf{\mathbf{u}} \bigr]
4203 * \underbrace{\bigl[ \mathsf{\mathbf{K}}_{uu} +
4204 * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
4205 * }_{\mathsf{\mathbf{K}}_{\textrm{con}}} d
4206 * \mathsf{\mathbf{u}}
4210 * \mathsf{\mathbf{F}}_{u}
4211 * - \mathsf{\mathbf{K}}_{u\widetilde{p}} \bigl[
4212 * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4213 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4215 * \overline{\mathsf{\mathbf{K}}}\mathsf{\mathbf{F}}_{\widetilde{p}}
4217 * \Bigr]}_{\mathsf{\mathbf{F}}_{\textrm{con}}}
4221 * \overline{\overline{\mathsf{\mathbf{K}}}} \dealcoloneq
4222 * \mathsf{\mathbf{K}}_{u\widetilde{p}}
4223 * \overline{\mathsf{\mathbf{K}}}
4224 * \mathsf{\mathbf{K}}_{\widetilde{p}u} \, .
4229 * At the top, we allocate two temporary vectors to help with the
4230 * static condensation, and variables to store the number of
4231 * linear solver iterations and the (hopefully converged) residual.
4234 * BlockVector<double> A(dofs_per_block);
4235 * BlockVector<double> B(dofs_per_block);
4240 * In the first step of this function, we solve for the incremental
4241 * displacement @f$d\mathbf{u}@f$. To this end, we perform static
4242 * condensation to make
4243 * @f$\mathsf{\mathbf{K}}_{\textrm{con}}
4244 * = \bigl[ \mathsf{\mathbf{K}}_{uu} +
4245 * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]@f$
4247 * @f$\mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
4248 * in the original @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
4249 * block. That is, we make @f$\mathsf{\mathbf{K}}_{\textrm{store}}@f$.
4258 * \mathsf{\mathbf{A}}_{\widetilde{J}}
4260 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4261 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4265 * tangent_matrix.block(p_dof, J_dof)
4266 * .vmult(A.block(J_dof), system_rhs.block(p_dof));
4270 * \mathsf{\mathbf{B}}_{\widetilde{J}}
4272 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4273 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4274 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4278 * tangent_matrix.block(J_dof, J_dof)
4279 * .vmult(B.block(J_dof), A.block(J_dof));
4283 * \mathsf{\mathbf{A}}_{\widetilde{J}}
4285 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4287 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4288 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4289 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4293 * A.block(J_dof) = system_rhs.block(J_dof);
4294 * A.block(J_dof) -= B.block(J_dof);
4298 * \mathsf{\mathbf{A}}_{\widetilde{J}}
4300 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4302 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4304 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4305 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4306 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4311 * tangent_matrix.block(p_dof, J_dof)
4312 * .Tvmult(A.block(p_dof), A.block(J_dof));
4316 * \mathsf{\mathbf{A}}_{u}
4318 * \mathsf{\mathbf{K}}_{u \widetilde{p}}
4319 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4321 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4323 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4324 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4325 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4330 * tangent_matrix.block(u_dof, p_dof)
4331 * .vmult(A.block(u_dof), A.block(p_dof));
4335 * \mathsf{\mathbf{F}}_{\text{con}}
4337 * \mathsf{\mathbf{F}}_{u}
4339 * \mathsf{\mathbf{K}}_{u \widetilde{p}}
4340 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4342 * \mathsf{\mathbf{F}}_{\widetilde{J}}
4344 * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4345 * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4346 * \mathsf{\mathbf{F}}_{\widetilde{p}}
4351 * system_rhs.block(u_dof) -= A.block(u_dof);
4353 * timer.enter_subsection("Linear solver");
4354 * std::cout << " SLV " << std::flush;
4355 * if (parameters.type_lin == "CG")
4357 * const auto solver_its = static_cast<unsigned int>(
4358 * tangent_matrix.block(u_dof, u_dof).m() *
4359 * parameters.max_iterations_lin);
4360 * const double tol_sol =
4361 * parameters.tol_lin * system_rhs.block(u_dof).l2_norm();
4363 * SolverControl solver_control(solver_its, tol_sol);
4365 * GrowingVectorMemory<Vector<double>> GVM;
4366 * SolverCG<Vector<double>> solver_CG(solver_control, GVM);
4377 *
preconditioner(parameters.preconditioner_type,
4378 *
parameters.preconditioner_relaxation);
4386 *
lin_it = solver_control.last_step();
4387 *
lin_res = solver_control.last_value();
4389 *
else if (parameters.type_lin ==
"Direct")
4407 *
Assert(
false, ExcMessage(
"Linear solver type not implemented"));
4409 *
timer.leave_subsection();
4420 *
timer.enter_subsection(
"Linear solver postprocessing");
4421 *
std::cout <<
" PP " << std::flush;
4459 *
A.block(
p_dof) *= -1.0;
4535 *
A.block(
J_dof) *= -1.0;
4577 *
timer.leave_subsection();
4581 *
std::cout <<
" ------ " << std::flush;
4583 *
timer.enter_subsection(
"Linear solver");
4584 *
std::cout <<
" SLV " << std::flush;
4586 *
if (parameters.type_lin ==
"CG")
4608 * blocks in
the RHS vector
4648 * diagonal), a Jacobi preconditioner
is suitable.
4657 *
parameters.max_iterations_lin),
4659 *
parameters.tol_lin);
4703 *
parameters.preconditioner_relaxation);
4708 *
parameters.max_iterations_lin),
4710 *
parameters.tol_lin);
4731 *
timer.leave_subsection();
4739 *
timer.enter_subsection(
"Linear solver postprocessing");
4740 *
std::cout <<
" PP " << std::flush;
4748 *
else if (parameters.type_lin ==
"Direct")
4767 *
std::cout <<
" -- " << std::flush;
4770 *
Assert(
false, ExcMessage(
"Linear solver type not implemented"));
4772 *
timer.leave_subsection();
4789 * <a name=
"step_44-Solidoutput_results"></a>
4805 *
template <
int dim>
4809 *
std::vector<DataComponentInterpretation::DataComponentInterpretation>
4810 *
data_component_interpretation(
4812 *
data_component_interpretation.push_back(
4814 *
data_component_interpretation.push_back(
4817 *
std::vector<std::string>
solution_name(dim,
"displacement");
4826 *
data_out.attach_dof_handler(dof_handler);
4830 *
data_component_interpretation);
4845 *
for (
unsigned int i = 0; i <
soln.size(); ++i)
4848 *
data_out.build_patches(
q_mapping, degree);
4850 *
std::ofstream output(
"solution-" + std::to_string(dim) +
"d-" +
4851 *
std::to_string(time.get_timestep()) +
".vtu");
4852 *
data_out.write_vtu(output);
4861 * <a name=
"step_44-Mainfunction"></a>
4869 *
using namespace Step44;
4873 *
const unsigned int dim = 3;
4877 *
catch (std::exception &exc)
4879 *
std::cerr << std::endl
4881 *
<<
"----------------------------------------------------"
4883 *
std::cerr <<
"Exception on processing: " << std::endl
4884 *
<< exc.what() << std::endl
4885 *
<<
"Aborting!" << std::endl
4886 *
<<
"----------------------------------------------------"
4893 *
std::cerr << std::endl
4895 *
<<
"----------------------------------------------------"
4897 *
std::cerr <<
"Unknown exception!" << std::endl
4898 *
<<
"Aborting!" << std::endl
4899 *
<<
"----------------------------------------------------"
4949 Reference volume: 1e-09
4969v /
V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
4988v /
V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
5004+---------------------------------------------+------------+------------+
5008+---------------------------------+-----------+------------+------------+
5011|
Linear solver | 43 | 9.248e+02s | 9.37e+01% |
5016+---------------------------------+-----------+------------+------------+
5048 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-displacement.png" alt=
"">
5054 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-pressure.png" alt=
"">
5060 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-dilatation.png" alt=
"">
5089 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-displacement.png" alt=
"">
5095 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-pressure.png" alt=
"">
5101 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-dilatation.png" alt=
"">
5116In terms of run-time, the @f$Q_2-DGPM_1-DGPM_1@f$ formulation tends to be more computationally expensive
5117than the @f$Q_1-DGPM_0-DGPM_0@f$ for a similar number of degrees-of-freedom
5118(produced by adding an extra grid refinement level for the lower-order interpolation).
5119This is shown in the graph below for a batch of tests run consecutively on a single 4-core (8-thread) machine.
5120The increase in computational time for the higher-order method is likely due to
5121the increased band-width required for the higher-order elements.
5122As previously mentioned, the use of a better solver and preconditioner may mitigate the
5123expense of using a higher-order formulation.
5124It was observed that for the given problem using the multithreaded Jacobi preconditioner can reduce the
5125computational runtime by up to 72% (for the worst case being a higher-order formulation with a large number
5126of degrees-of-freedom) in comparison to the single-thread SSOR preconditioner.
5134 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.Normalised_runtime.png" alt=
"">
5152 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.2d-gr_2.png" alt=
"">
5158 <
img src=
"https://www.dealii.org/images/steps/developer/step-44.2d-gr_5.png" alt=
"">
5166<a name=
"step-44-extensions"></a>
5167<a name=
"step_44-Possibilitiesforextensions"></a><
h3>Possibilities
for extensions</
h3>
5218<a name=
"step_44-PlainProg"></a>
numbers::NumberTraits< Number >::real_type norm() const
std::conditional_t< rank_==1, std::array< Number, dim >, std::array< Tensor< rank_ - 1, dim, Number >, dim > > values
#define Assert(cond, exc)
#define AssertDimension(dim1, dim2)
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
LinearOperator< Range_2, Domain_2, Payload > schur_complement(const LinearOperator< Domain_1, Range_1, Payload > &A_inv, const LinearOperator< Range_1, Domain_2, Payload > &B, const LinearOperator< Range_2, Domain_1, Payload > &C, const LinearOperator< Range_2, Domain_2, Payload > &D)
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
void loop(IteratorType begin, std_cxx20::type_identity_t< IteratorType > end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(std_cxx20::type_identity_t< DOFINFO > &, std_cxx20::type_identity_t< DOFINFO > &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, AssemblerType &assembler, const LoopControl &lctrl=LoopControl())
Task< RT > new_task(const std::function< RT()> &function)
#define DEAL_II_ASSERT_UNREACHABLE()
std::vector< index_type > data
std::vector< value_type > split(const typename ::Triangulation< dim, spacedim >::cell_iterator &parent, const value_type parent_value)
CGAL::Exact_predicates_exact_constructions_kernel_with_sqrt K
@ component_is_part_of_vector
Expression sign(const Expression &x)
@ matrix
Contents is actually a matrix.
@ symmetric
Matrix is symmetric.
@ diagonal
Matrix is diagonal.
@ general
No special properties.
constexpr types::blas_int zero
constexpr types::blas_int one
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.)
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
VectorType::value_type * end(VectorType &V)
bool check(const ConstraintKinds kind_in, const unsigned int dim)
DEAL_II_HOST constexpr TableIndices< 2 > merge(const TableIndices< 2 > &previous_indices, const unsigned int new_index, const unsigned int position)
void copy(const T *begin, const T *end, U *dest)
int(&) functions(const void *v1, const void *v2)
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
const ::parallel::distributed::Triangulation< dim, spacedim > * triangulation