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-44.h
Go to the documentation of this file.
1 )
1898  * {}
1899  *
1900  * @endcode
1901  *
1902  * Solve for the displacement using a Newton-Raphson method. We break this
1903  * function into the nonlinear loop and the function that solves the
1904  * linearized Newton-Raphson step:
1905  *
1906  * @code
1907  * void solve_nonlinear_timestep(BlockVector<double> &solution_delta);
1908  *
1909  * std::pair<unsigned int, double>
1910  * solve_linear_system(BlockVector<double> &newton_update);
1911  *
1912  * @endcode
1913  *
1914  * Solution retrieval as well as post-processing and writing data to file :
1915  *
1916  * @code
1918  * get_total_solution(const BlockVector<double> &solution_delta) const;
1919  *
1920  * void output_results() const;
1921  *
1922  * @endcode
1923  *
1924  * Finally, some member variables that describe the current state: A
1925  * collection of the parameters used to describe the problem setup...
1926  *
1927  * @code
1928  * Parameters::AllParameters parameters;
1929  *
1930  * @endcode
1931  *
1932  * ...the volume of the reference configuration...
1933  *
1934  * @code
1935  * double vol_reference;
1936  *
1937  * @endcode
1938  *
1939  * ...and description of the geometry on which the problem is solved:
1940  *
1941  * @code
1943  *
1944  * @endcode
1945  *
1946  * Also, keep track of the current time and the time spent evaluating
1947  * certain functions
1948  *
1949  * @code
1950  * Time time;
1951  * mutable TimerOutput timer;
1952  *
1953  * @endcode
1954  *
1955  * A storage object for quadrature point information. As opposed to
1956  * @ref step_18 "step-18", deal.II's native quadrature point data manager is employed
1957  * here.
1958  *
1959  * @code
1960  * CellDataStorage<typename Triangulation<dim>::cell_iterator,
1961  * PointHistory<dim>>
1962  * quadrature_point_history;
1963  *
1964  * @endcode
1965  *
1966  * A description of the finite-element system including the displacement
1967  * polynomial degree, the degree-of-freedom handler, number of DoFs per
1968  * cell and the extractor objects used to retrieve information from the
1969  * solution vectors:
1970  *
1971  * @code
1972  * const unsigned int degree;
1973  * const FESystem<dim> fe;
1974  * DoFHandler<dim> dof_handler;
1975  * const unsigned int dofs_per_cell;
1976  * const FEValuesExtractors::Vector u_fe;
1977  * const FEValuesExtractors::Scalar p_fe;
1978  * const FEValuesExtractors::Scalar J_fe;
1979  *
1980  * @endcode
1981  *
1982  * Description of how the block-system is arranged. There are 3 blocks,
1983  * the first contains a vector DOF @f$\mathbf{u}@f$ while the other two
1984  * describe scalar DOFs, @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$.
1985  *
1986  * @code
1987  * static const unsigned int n_blocks = 3;
1988  * static const unsigned int n_components = dim + 2;
1989  * static const unsigned int first_u_component = 0;
1990  * static const unsigned int p_component = dim;
1991  * static const unsigned int J_component = dim + 1;
1992  *
1993  * enum
1994  * {
1995  * u_dof = 0,
1996  * p_dof = 1,
1997  * J_dof = 2
1998  * };
1999  *
2000  * std::vector<types::global_dof_index> dofs_per_block;
2001  * std::vector<types::global_dof_index> element_indices_u;
2002  * std::vector<types::global_dof_index> element_indices_p;
2003  * std::vector<types::global_dof_index> element_indices_J;
2004  *
2005  * @endcode
2006  *
2007  * Rules for Gauss-quadrature on both the cell and faces. The number of
2008  * quadrature points on both cells and faces is recorded.
2009  *
2010  * @code
2011  * const QGauss<dim> qf_cell;
2012  * const QGauss<dim - 1> qf_face;
2013  * const unsigned int n_q_points;
2014  * const unsigned int n_q_points_f;
2015  *
2016  * @endcode
2017  *
2018  * Objects that store the converged solution and right-hand side vectors,
2019  * as well as the tangent matrix. There is an AffineConstraints object used
2020  * to keep track of constraints. We make use of a sparsity pattern
2021  * designed for a block system.
2022  *
2023  * @code
2024  * AffineConstraints<double> constraints;
2025  * BlockSparsityPattern sparsity_pattern;
2026  * BlockSparseMatrix<double> tangent_matrix;
2027  * BlockVector<double> system_rhs;
2028  * BlockVector<double> solution_n;
2029  *
2030  * @endcode
2031  *
2032  * Then define a number of variables to store norms and update norms and
2033  * normalization factors.
2034  *
2035  * @code
2036  * struct Errors
2037  * {
2038  * Errors()
2039  * : norm(1.0)
2040  * , u(1.0)
2041  * , p(1.0)
2042  * , J(1.0)
2043  * {}
2044  *
2045  * void reset()
2046  * {
2047  * norm = 1.0;
2048  * u = 1.0;
2049  * p = 1.0;
2050  * J = 1.0;
2051  * }
2052  * void normalize(const Errors &rhs)
2053  * {
2054  * if (rhs.norm != 0.0)
2055  * norm /= rhs.norm;
2056  * if (rhs.u != 0.0)
2057  * u /= rhs.u;
2058  * if (rhs.p != 0.0)
2059  * p /= rhs.p;
2060  * if (rhs.J != 0.0)
2061  * J /= rhs.J;
2062  * }
2063  *
2064  * double norm, u, p, J;
2065  * };
2066  *
2067  * Errors error_residual, error_residual_0, error_residual_norm, error_update,
2068  * error_update_0, error_update_norm;
2069  *
2070  * @endcode
2071  *
2072  * Methods to calculate error measures
2073  *
2074  * @code
2075  * void get_error_residual(Errors &error_residual);
2076  *
2077  * void get_error_update(const BlockVector<double> &newton_update,
2078  * Errors & error_update);
2079  *
2080  * std::pair<double, double> get_error_dilation() const;
2081  *
2082  * @endcode
2083  *
2084  * Compute the volume in the spatial configuration
2085  *
2086  * @code
2087  * double compute_vol_current() const;
2088  *
2089  * @endcode
2090  *
2091  * Print information to screen in a pleasing way...
2092  *
2093  * @code
2094  * static void print_conv_header();
2095  *
2096  * void print_conv_footer();
2097  * };
2098  *
2099  * @endcode
2100  *
2101  *
2102  * <a name="ImplementationofthecodeSolidcodeclass"></a>
2103  * <h3>Implementation of the <code>Solid</code> class</h3>
2104  *
2105 
2106  *
2107  *
2108  * <a name="Publicinterface"></a>
2109  * <h4>Public interface</h4>
2110  *
2111 
2112  *
2113  * We initialize the Solid class using data extracted from the parameter file.
2114  *
2115  * @code
2116  * template <int dim>
2117  * Solid<dim>::Solid(const std::string &input_file)
2118  * : parameters(input_file)
2119  * , vol_reference(0.)
2120  * , triangulation(Triangulation<dim>::maximum_smoothing)
2121  * , time(parameters.end_time, parameters.delta_t)
2122  * , timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
2123  * , degree(parameters.poly_degree)
2124  * ,
2125  * @endcode
2126  *
2127  * The Finite Element System is composed of dim continuous displacement
2128  * DOFs, and discontinuous pressure and dilatation DOFs. In an attempt to
2129  * satisfy the Babuska-Brezzi or LBB stability conditions (see Hughes
2130  * (2000)), we setup a @f$Q_n \times DGPM_{n-1} \times DGPM_{n-1}@f$
2131  * system. @f$Q_2 \times DGPM_1 \times DGPM_1@f$ elements satisfy this
2132  * condition, while @f$Q_1 \times DGPM_0 \times DGPM_0@f$ elements do
2133  * not. However, it has been shown that the latter demonstrate good
2134  * convergence characteristics nonetheless.
2135  *
2136  * @code
2137  * fe(FE_Q<dim>(parameters.poly_degree),
2138  * dim, // displacement
2139  * FE_DGPMonomial<dim>(parameters.poly_degree - 1),
2140  * 1, // pressure
2141  * FE_DGPMonomial<dim>(parameters.poly_degree - 1),
2142  * 1)
2143  * , // dilatation
2144  * dof_handler(triangulation)
2145  * , dofs_per_cell(fe.dofs_per_cell)
2146  * , u_fe(first_u_component)
2147  * , p_fe(p_component)
2148  * , J_fe(J_component)
2149  * , dofs_per_block(n_blocks)
2150  * , qf_cell(parameters.quad_order)
2151  * , qf_face(parameters.quad_order)
2152  * , n_q_points(qf_cell.size())
2153  * , n_q_points_f(qf_face.size())
2154  * {
2155  * Assert(dim == 2 || dim == 3,
2156  * ExcMessage("This problem only works in 2 or 3 space dimensions."));
2157  * determine_component_extractors();
2158  * }
2159  *
2160  *
2161  * @endcode
2162  *
2163  * In solving the quasi-static problem, the time becomes a loading parameter,
2164  * i.e. we increasing the loading linearly with time, making the two concepts
2165  * interchangeable. We choose to increment time linearly using a constant time
2166  * step size.
2167  *
2168 
2169  *
2170  * We start the function with preprocessing, setting the initial dilatation
2171  * values, and then output the initial grid before starting the simulation
2172  * proper with the first time (and loading)
2173  * increment.
2174  *
2175 
2176  *
2177  * Care must be taken (or at least some thought given) when imposing the
2178  * constraint @f$\widetilde{J}=1@f$ on the initial solution field. The constraint
2179  * corresponds to the determinant of the deformation gradient in the
2180  * undeformed configuration, which is the identity tensor. We use
2181  * FE_DGPMonomial bases to interpolate the dilatation field, thus we can't
2182  * simply set the corresponding dof to unity as they correspond to the
2183  * monomial coefficients. Thus we use the VectorTools::project function to do
2184  * the work for us. The VectorTools::project function requires an argument
2185  * indicating the hanging node constraints. We have none in this program
2186  * So we have to create a constraint object. In its original state, constraint
2187  * objects are unsorted, and have to be sorted (using the
2188  * AffineConstraints::close function) before they can be used. Have a look at
2189  * @ref step_21 "step-21" for more information. We only need to enforce the initial condition
2190  * on the dilatation. In order to do this, we make use of a
2191  * ComponentSelectFunction which acts as a mask and sets the J_component of
2192  * n_components to 1. This is exactly what we want. Have a look at its usage
2193  * in @ref step_20 "step-20" for more information.
2194  *
2195  * @code
2196  * template <int dim>
2197  * void Solid<dim>::run()
2198  * {
2199  * make_grid();
2200  * system_setup();
2201  * {
2202  * AffineConstraints<double> constraints;
2203  * constraints.close();
2204  *
2205  * const ComponentSelectFunction<dim> J_mask(J_component, n_components);
2206  *
2208  * dof_handler, constraints, QGauss<dim>(degree + 2), J_mask, solution_n);
2209  * }
2210  * output_results();
2211  * time.increment();
2212  *
2213  * @endcode
2214  *
2215  * We then declare the incremental solution update @f$\varDelta
2216  * \mathbf{\Xi} \dealcoloneq \{\varDelta \mathbf{u},\varDelta \widetilde{p},
2217  * \varDelta \widetilde{J} \}@f$ and start the loop over the time domain.
2218  *
2219 
2220  *
2221  * At the beginning, we reset the solution update for this time step...
2222  *
2223  * @code
2224  * BlockVector<double> solution_delta(dofs_per_block);
2225  * while (time.current() < time.end())
2226  * {
2227  * solution_delta = 0.0;
2228  *
2229  * @endcode
2230  *
2231  * ...solve the current time step and update total solution vector
2232  * @f$\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
2233  * \varDelta \mathbf{\Xi}@f$...
2234  *
2235  * @code
2236  * solve_nonlinear_timestep(solution_delta);
2237  * solution_n += solution_delta;
2238  *
2239  * @endcode
2240  *
2241  * ...and plot the results before moving on happily to the next time
2242  * step:
2243  *
2244  * @code
2245  * output_results();
2246  * time.increment();
2247  * }
2248  * }
2249  *
2250  *
2251  * @endcode
2252  *
2253  *
2254  * <a name="Privateinterface"></a>
2255  * <h3>Private interface</h3>
2256  *
2257 
2258  *
2259  *
2260  * <a name="Threadingbuildingblocksstructures"></a>
2261  * <h4>Threading-building-blocks structures</h4>
2262  *
2263 
2264  *
2265  * The first group of private member functions is related to parallelization.
2266  * We use the Threading Building Blocks library (TBB) to perform as many
2267  * computationally intensive distributed tasks as possible. In particular, we
2268  * assemble the tangent matrix and right hand side vector, the static
2269  * condensation contributions, and update data stored at the quadrature points
2270  * using TBB. Our main tool for this is the WorkStream class (see the @ref
2271  * threads module for more information).
2272  *
2273 
2274  *
2275  * Firstly we deal with the tangent matrix assembly structures. The
2276  * PerTaskData object stores local contributions.
2277  *
2278  * @code
2279  * template <int dim>
2280  * struct Solid<dim>::PerTaskData_K
2281  * {
2283  * std::vector<types::global_dof_index> local_dof_indices;
2284  *
2285  * PerTaskData_K(const unsigned int dofs_per_cell)
2286  * : cell_matrix(dofs_per_cell, dofs_per_cell)
2287  * , local_dof_indices(dofs_per_cell)
2288  * {}
2289  *
2290  * void reset()
2291  * {
2292  * cell_matrix = 0.0;
2293  * }
2294  * };
2295  *
2296  *
2297  * @endcode
2298  *
2299  * On the other hand, the ScratchData object stores the larger objects such as
2300  * the shape-function values array (<code>Nx</code>) and a shape function
2301  * gradient and symmetric gradient vector which we will use during the
2302  * assembly.
2303  *
2304  * @code
2305  * template <int dim>
2306  * struct Solid<dim>::ScratchData_K
2307  * {
2308  * FEValues<dim> fe_values;
2309  *
2310  * std::vector<std::vector<double>> Nx;
2311  * std::vector<std::vector<Tensor<2, dim>>> grad_Nx;
2312  * std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
2313  *
2314  * ScratchData_K(const FiniteElement<dim> &fe_cell,
2315  * const QGauss<dim> & qf_cell,
2316  * const UpdateFlags uf_cell)
2317  * : fe_values(fe_cell, qf_cell, uf_cell)
2318  * , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
2319  * , grad_Nx(qf_cell.size(),
2320  * std::vector<Tensor<2, dim>>(fe_cell.dofs_per_cell))
2321  * , symm_grad_Nx(qf_cell.size(),
2322  * std::vector<SymmetricTensor<2, dim>>(
2323  * fe_cell.dofs_per_cell))
2324  * {}
2325  *
2326  * ScratchData_K(const ScratchData_K &rhs)
2327  * : fe_values(rhs.fe_values.get_fe(),
2328  * rhs.fe_values.get_quadrature(),
2329  * rhs.fe_values.get_update_flags())
2330  * , Nx(rhs.Nx)
2331  * , grad_Nx(rhs.grad_Nx)
2332  * , symm_grad_Nx(rhs.symm_grad_Nx)
2333  * {}
2334  *
2335  * void reset()
2336  * {
2337  * const unsigned int n_q_points = Nx.size();
2338  * const unsigned int n_dofs_per_cell = Nx[0].size();
2339  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2340  * {
2341  * Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
2342  * Assert(grad_Nx[q_point].size() == n_dofs_per_cell,
2343  * ExcInternalError());
2344  * Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
2345  * ExcInternalError());
2346  * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2347  * {
2348  * Nx[q_point][k] = 0.0;
2349  * grad_Nx[q_point][k] = 0.0;
2350  * symm_grad_Nx[q_point][k] = 0.0;
2351  * }
2352  * }
2353  * }
2354  * };
2355  *
2356  * @endcode
2357  *
2358  * Next, the same approach is used for the right-hand side assembly. The
2359  * PerTaskData object again stores local contributions and the ScratchData
2360  * object the shape function object and precomputed values vector:
2361  *
2362  * @code
2363  * template <int dim>
2364  * struct Solid<dim>::PerTaskData_RHS
2365  * {
2366  * Vector<double> cell_rhs;
2367  * std::vector<types::global_dof_index> local_dof_indices;
2368  *
2369  * PerTaskData_RHS(const unsigned int dofs_per_cell)
2370  * : cell_rhs(dofs_per_cell)
2371  * , local_dof_indices(dofs_per_cell)
2372  * {}
2373  *
2374  * void reset()
2375  * {
2376  * cell_rhs = 0.0;
2377  * }
2378  * };
2379  *
2380  *
2381  * template <int dim>
2382  * struct Solid<dim>::ScratchData_RHS
2383  * {
2384  * FEValues<dim> fe_values;
2385  * FEFaceValues<dim> fe_face_values;
2386  *
2387  * std::vector<std::vector<double>> Nx;
2388  * std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
2389  *
2390  * ScratchData_RHS(const FiniteElement<dim> &fe_cell,
2391  * const QGauss<dim> & qf_cell,
2392  * const UpdateFlags uf_cell,
2393  * const QGauss<dim - 1> & qf_face,
2394  * const UpdateFlags uf_face)
2395  * : fe_values(fe_cell, qf_cell, uf_cell)
2396  * , fe_face_values(fe_cell, qf_face, uf_face)
2397  * , Nx(qf_cell.size(), std::vector<double>(fe_cell.dofs_per_cell))
2398  * , symm_grad_Nx(qf_cell.size(),
2399  * std::vector<SymmetricTensor<2, dim>>(
2400  * fe_cell.dofs_per_cell))
2401  * {}
2402  *
2403  * ScratchData_RHS(const ScratchData_RHS &rhs)
2404  * : fe_values(rhs.fe_values.get_fe(),
2405  * rhs.fe_values.get_quadrature(),
2406  * rhs.fe_values.get_update_flags())
2407  * , fe_face_values(rhs.fe_face_values.get_fe(),
2408  * rhs.fe_face_values.get_quadrature(),
2409  * rhs.fe_face_values.get_update_flags())
2410  * , Nx(rhs.Nx)
2411  * , symm_grad_Nx(rhs.symm_grad_Nx)
2412  * {}
2413  *
2414  * void reset()
2415  * {
2416  * const unsigned int n_q_points = Nx.size();
2417  * const unsigned int n_dofs_per_cell = Nx[0].size();
2418  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2419  * {
2420  * Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
2421  * Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
2422  * ExcInternalError());
2423  * for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
2424  * {
2425  * Nx[q_point][k] = 0.0;
2426  * symm_grad_Nx[q_point][k] = 0.0;
2427  * }
2428  * }
2429  * }
2430  * };
2431  *
2432  * @endcode
2433  *
2434  * Then we define structures to assemble the statically condensed tangent
2435  * matrix. Recall that we wish to solve for a displacement-based formulation.
2436  * We do the condensation at the element level as the @f$\widetilde{p}@f$ and
2437  * @f$\widetilde{J}@f$ fields are element-wise discontinuous. As these operations
2438  * are matrix-based, we need to setup a number of matrices to store the local
2439  * contributions from a number of the tangent matrix sub-blocks. We place
2440  * these in the PerTaskData struct.
2441  *
2442 
2443  *
2444  * We choose not to reset any data in the <code>reset()</code> function as the
2445  * matrix extraction and replacement tools will take care of this
2446  *
2447  * @code
2448  * template <int dim>
2449  * struct Solid<dim>::PerTaskData_SC
2450  * {
2452  * std::vector<types::global_dof_index> local_dof_indices;
2453  *
2454  * FullMatrix<double> k_orig;
2455  * FullMatrix<double> k_pu;
2456  * FullMatrix<double> k_pJ;
2457  * FullMatrix<double> k_JJ;
2458  * FullMatrix<double> k_pJ_inv;
2459  * FullMatrix<double> k_bbar;
2461  * FullMatrix<double> B;
2463  *
2464  * PerTaskData_SC(const unsigned int dofs_per_cell,
2465  * const unsigned int n_u,
2466  * const unsigned int n_p,
2467  * const unsigned int n_J)
2468  * : cell_matrix(dofs_per_cell, dofs_per_cell)
2469  * , local_dof_indices(dofs_per_cell)
2470  * , k_orig(dofs_per_cell, dofs_per_cell)
2471  * , k_pu(n_p, n_u)
2472  * , k_pJ(n_p, n_J)
2473  * , k_JJ(n_J, n_J)
2474  * , k_pJ_inv(n_p, n_J)
2475  * , k_bbar(n_u, n_u)
2476  * , A(n_J, n_u)
2477  * , B(n_J, n_u)
2478  * , C(n_p, n_u)
2479  * {}
2480  *
2481  * void reset()
2482  * {}
2483  * };
2484  *
2485  *
2486  * @endcode
2487  *
2488  * The ScratchData object for the operations we wish to perform here is empty
2489  * since we need no temporary data, but it still needs to be defined for the
2490  * current implementation of TBB in deal.II. So we create a dummy struct for
2491  * this purpose.
2492  *
2493  * @code
2494  * template <int dim>
2495  * struct Solid<dim>::ScratchData_SC
2496  * {
2497  * void reset()
2498  * {}
2499  * };
2500  *
2501  *
2502  * @endcode
2503  *
2504  * And finally we define the structures to assist with updating the quadrature
2505  * point information. Similar to the SC assembly process, we do not need the
2506  * PerTaskData object (since there is nothing to store here) but must define
2507  * one nonetheless. Note that this is because for the operation that we have
2508  * here -- updating the data on quadrature points -- the operation is purely
2509  * local: the things we do on every cell get consumed on every cell, without
2510  * any global aggregation operation as is usually the case when using the
2511  * WorkStream class. The fact that we still have to define a per-task data
2512  * structure points to the fact that the WorkStream class may be ill-suited to
2513  * this operation (we could, in principle simply create a new task using
2514  * Threads::new_task for each cell) but there is not much harm done to doing
2515  * it this way anyway.
2516  * Furthermore, should there be different material models associated with a
2517  * quadrature point, requiring varying levels of computational expense, then
2518  * the method used here could be advantageous.
2519  *
2520  * @code
2521  * template <int dim>
2522  * struct Solid<dim>::PerTaskData_UQPH
2523  * {
2524  * void reset()
2525  * {}
2526  * };
2527  *
2528  *
2529  * @endcode
2530  *
2531  * The ScratchData object will be used to store an alias for the solution
2532  * vector so that we don't have to copy this large data structure. We then
2533  * define a number of vectors to extract the solution values and gradients at
2534  * the quadrature points.
2535  *
2536  * @code
2537  * template <int dim>
2538  * struct Solid<dim>::ScratchData_UQPH
2539  * {
2540  * const BlockVector<double> &solution_total;
2541  *
2542  * std::vector<Tensor<2, dim>> solution_grads_u_total;
2543  * std::vector<double> solution_values_p_total;
2544  * std::vector<double> solution_values_J_total;
2545  *
2546  * FEValues<dim> fe_values;
2547  *
2548  * ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
2549  * const QGauss<dim> & qf_cell,
2550  * const UpdateFlags uf_cell,
2551  * const BlockVector<double> &solution_total)
2552  * : solution_total(solution_total)
2553  * , solution_grads_u_total(qf_cell.size())
2554  * , solution_values_p_total(qf_cell.size())
2555  * , solution_values_J_total(qf_cell.size())
2556  * , fe_values(fe_cell, qf_cell, uf_cell)
2557  * {}
2558  *
2559  * ScratchData_UQPH(const ScratchData_UQPH &rhs)
2560  * : solution_total(rhs.solution_total)
2561  * , solution_grads_u_total(rhs.solution_grads_u_total)
2562  * , solution_values_p_total(rhs.solution_values_p_total)
2563  * , solution_values_J_total(rhs.solution_values_J_total)
2564  * , fe_values(rhs.fe_values.get_fe(),
2565  * rhs.fe_values.get_quadrature(),
2566  * rhs.fe_values.get_update_flags())
2567  * {}
2568  *
2569  * void reset()
2570  * {
2571  * const unsigned int n_q_points = solution_grads_u_total.size();
2572  * for (unsigned int q = 0; q < n_q_points; ++q)
2573  * {
2574  * solution_grads_u_total[q] = 0.0;
2575  * solution_values_p_total[q] = 0.0;
2576  * solution_values_J_total[q] = 0.0;
2577  * }
2578  * }
2579  * };
2580  *
2581  *
2582  * @endcode
2583  *
2584  *
2585  * <a name="Solidmake_grid"></a>
2586  * <h4>Solid::make_grid</h4>
2587  *
2588 
2589  *
2590  * On to the first of the private member functions. Here we create the
2591  * triangulation of the domain, for which we choose the scaled cube with each
2592  * face given a boundary ID number. The grid must be refined at least once
2593  * for the indentation problem.
2594  *
2595 
2596  *
2597  * We then determine the volume of the reference configuration and print it
2598  * for comparison:
2599  *
2600  * @code
2601  * template <int dim>
2602  * void Solid<dim>::make_grid()
2603  * {
2604  * GridGenerator::hyper_rectangle(
2605  * triangulation,
2606  * (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
2607  * (dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
2608  * true);
2609  * GridTools::scale(parameters.scale, triangulation);
2610  * triangulation.refine_global(std::max(1U, parameters.global_refinement));
2611  *
2612  * vol_reference = GridTools::volume(triangulation);
2613  * std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
2614  *
2615  * @endcode
2616  *
2617  * Since we wish to apply a Neumann BC to a patch on the top surface, we
2618  * must find the cell faces in this part of the domain and mark them with
2619  * a distinct boundary ID number. The faces we are looking for are on the
2620  * +y surface and will get boundary ID 6 (zero through five are already
2621  * used when creating the six faces of the cube domain):
2622  *
2623  * @code
2624  * for (const auto &cell : triangulation.active_cell_iterators())
2625  * for (const auto &face : cell->face_iterators())
2626  * {
2627  * if (face->at_boundary() == true &&
2628  * face->center()[1] == 1.0 * parameters.scale)
2629  * {
2630  * if (dim == 3)
2631  * {
2632  * if (face->center()[0] < 0.5 * parameters.scale &&
2633  * face->center()[2] < 0.5 * parameters.scale)
2634  * face->set_boundary_id(6);
2635  * }
2636  * else
2637  * {
2638  * if (face->center()[0] < 0.5 * parameters.scale)
2639  * face->set_boundary_id(6);
2640  * }
2641  * }
2642  * }
2643  * }
2644  *
2645  *
2646  * @endcode
2647  *
2648  *
2649  * <a name="Solidsystem_setup"></a>
2650  * <h4>Solid::system_setup</h4>
2651  *
2652 
2653  *
2654  * Next we describe how the FE system is setup. We first determine the number
2655  * of components per block. Since the displacement is a vector component, the
2656  * first dim components belong to it, while the next two describe scalar
2657  * pressure and dilatation DOFs.
2658  *
2659  * @code
2660  * template <int dim>
2661  * void Solid<dim>::system_setup()
2662  * {
2663  * timer.enter_subsection("Setup system");
2664  *
2665  * std::vector<unsigned int> block_component(n_components,
2666  * u_dof); // Displacement
2667  * block_component[p_component] = p_dof; // Pressure
2668  * block_component[J_component] = J_dof; // Dilatation
2669  *
2670  * @endcode
2671  *
2672  * The DOF handler is then initialized and we renumber the grid in an
2673  * efficient manner. We also record the number of DOFs per block.
2674  *
2675  * @code
2676  * dof_handler.distribute_dofs(fe);
2677  * DoFRenumbering::Cuthill_McKee(dof_handler);
2678  * DoFRenumbering::component_wise(dof_handler, block_component);
2679  *
2680  * dofs_per_block =
2681  * DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
2682  *
2683  * std::cout << "Triangulation:"
2684  * << "\n\t Number of active cells: "
2685  * << triangulation.n_active_cells()
2686  * << "\n\t Number of degrees of freedom: " << dof_handler.n_dofs()
2687  * << std::endl;
2688  *
2689  * @endcode
2690  *
2691  * Setup the sparsity pattern and tangent matrix
2692  *
2693  * @code
2694  * tangent_matrix.clear();
2695  * {
2696  * const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
2697  * const types::global_dof_index n_dofs_p = dofs_per_block[p_dof];
2698  * const types::global_dof_index n_dofs_J = dofs_per_block[J_dof];
2699  *
2700  * BlockDynamicSparsityPattern dsp(n_blocks, n_blocks);
2701  *
2702  * dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
2703  * dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
2704  * dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
2705  *
2706  * dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
2707  * dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
2708  * dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
2709  *
2710  * dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
2711  * dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
2712  * dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
2713  * dsp.collect_sizes();
2714  *
2715  * @endcode
2716  *
2717  * The global system matrix initially has the following structure
2718  * @f{align*}
2719  * \underbrace{\begin{bmatrix}
2720  * \mathsf{\mathbf{K}}_{uu} & \mathsf{\mathbf{K}}_{u\widetilde{p}} &
2721  * \mathbf{0}
2722  * \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
2723  * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}
2724  * \\ \mathbf{0} & \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
2725  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
2726  * \end{bmatrix}}_{\mathsf{\mathbf{K}}(\mathbf{\Xi}_{\textrm{i}})}
2727  * \underbrace{\begin{bmatrix}
2728  * d \mathsf{u}
2729  * \\ d \widetilde{\mathsf{\mathbf{p}}}
2730  * \\ d \widetilde{\mathsf{\mathbf{J}}}
2731  * \end{bmatrix}}_{d \mathbf{\Xi}}
2732  * =
2733  * \underbrace{\begin{bmatrix}
2734  * \mathsf{\mathbf{F}}_{u}(\mathbf{u}_{\textrm{i}})
2735  * \\ \mathsf{\mathbf{F}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}})
2736  * \\ \mathsf{\mathbf{F}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}})
2737  * \end{bmatrix}}_{ \mathsf{\mathbf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, .
2738  * @f}
2739  * We optimize the sparsity pattern to reflect this structure
2740  * and prevent unnecessary data creation for the right-diagonal
2741  * block components.
2742  *
2743  * @code
2744  * Table<2, DoFTools::Coupling> coupling(n_components, n_components);
2745  * for (unsigned int ii = 0; ii < n_components; ++ii)
2746  * for (unsigned int jj = 0; jj < n_components; ++jj)
2747  * if (((ii < p_component) && (jj == J_component)) ||
2748  * ((ii == J_component) && (jj < p_component)) ||
2749  * ((ii == p_component) && (jj == p_component)))
2750  * coupling[ii][jj] = DoFTools::none;
2751  * else
2752  * coupling[ii][jj] = DoFTools::always;
2753  * DoFTools::make_sparsity_pattern(
2754  * dof_handler, coupling, dsp, constraints, false);
2755  * sparsity_pattern.copy_from(dsp);
2756  * }
2757  *
2758  * tangent_matrix.reinit(sparsity_pattern);
2759  *
2760  * @endcode
2761  *
2762  * We then set up storage vectors
2763  *
2764  * @code
2765  * system_rhs.reinit(dofs_per_block);
2766  * system_rhs.collect_sizes();
2767  *
2768  * solution_n.reinit(dofs_per_block);
2769  * solution_n.collect_sizes();
2770  *
2771  * @endcode
2772  *
2773  * ...and finally set up the quadrature
2774  * point history:
2775  *
2776  * @code
2777  * setup_qph();
2778  *
2779  * timer.leave_subsection();
2780  * }
2781  *
2782  *
2783  * @endcode
2784  *
2785  *
2786  * <a name="Soliddetermine_component_extractors"></a>
2787  * <h4>Solid::determine_component_extractors</h4>
2788  * Next we compute some information from the FE system that describes which
2789  * local element DOFs are attached to which block component. This is used
2790  * later to extract sub-blocks from the global matrix.
2791  *
2792 
2793  *
2794  * In essence, all we need is for the FESystem object to indicate to which
2795  * block component a DOF on the reference cell is attached to. Currently, the
2796  * interpolation fields are setup such that 0 indicates a displacement DOF, 1
2797  * a pressure DOF and 2 a dilatation DOF.
2798  *
2799  * @code
2800  * template <int dim>
2801  * void Solid<dim>::determine_component_extractors()
2802  * {
2803  * element_indices_u.clear();
2804  * element_indices_p.clear();
2805  * element_indices_J.clear();
2806  *
2807  * for (unsigned int k = 0; k < fe.dofs_per_cell; ++k)
2808  * {
2809  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
2810  * if (k_group == u_dof)
2811  * element_indices_u.push_back(k);
2812  * else if (k_group == p_dof)
2813  * element_indices_p.push_back(k);
2814  * else if (k_group == J_dof)
2815  * element_indices_J.push_back(k);
2816  * else
2817  * {
2818  * Assert(k_group <= J_dof, ExcInternalError());
2819  * }
2820  * }
2821  * }
2822  *
2823  * @endcode
2824  *
2825  *
2826  * <a name="Solidsetup_qph"></a>
2827  * <h4>Solid::setup_qph</h4>
2828  * The method used to store quadrature information is already described in
2829  * @ref step_18 "step-18". Here we implement a similar setup for a SMP machine.
2830  *
2831 
2832  *
2833  * Firstly the actual QPH data objects are created. This must be done only
2834  * once the grid is refined to its finest level.
2835  *
2836  * @code
2837  * template <int dim>
2838  * void Solid<dim>::setup_qph()
2839  * {
2840  * std::cout << " Setting up quadrature point data..." << std::endl;
2841  *
2842  * quadrature_point_history.initialize(triangulation.begin_active(),
2843  * triangulation.end(),
2844  * n_q_points);
2845  *
2846  * @endcode
2847  *
2848  * Next we setup the initial quadrature point data.
2849  * Note that when the quadrature point data is retrieved,
2850  * it is returned as a vector of smart pointers.
2851  *
2852  * @code
2853  * for (const auto &cell : triangulation.active_cell_iterators())
2854  * {
2855  * const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
2856  * quadrature_point_history.get_data(cell);
2857  * Assert(lqph.size() == n_q_points, ExcInternalError());
2858  *
2859  * for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
2860  * lqph[q_point]->setup_lqp(parameters);
2861  * }
2862  * }
2863  *
2864  * @endcode
2865  *
2866  *
2867  * <a name="Solidupdate_qph_incremental"></a>
2868  * <h4>Solid::update_qph_incremental</h4>
2869  * As the update of QP information occurs frequently and involves a number of
2870  * expensive operations, we define a multithreaded approach to distributing
2871  * the task across a number of CPU cores.
2872  *
2873 
2874  *
2875  * To start this, we first we need to obtain the total solution as it stands
2876  * at this Newton increment and then create the initial copy of the scratch
2877  * and copy data objects:
2878  *
2879  * @code
2880  * template <int dim>
2881  * void
2882  * Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
2883  * {
2884  * timer.enter_subsection("Update QPH data");
2885  * std::cout << " UQPH " << std::flush;
2886  *
2887  * const BlockVector<double> solution_total(
2888  * get_total_solution(solution_delta));
2889  *
2890  * const UpdateFlags uf_UQPH(update_values | update_gradients);
2891  * PerTaskData_UQPH per_task_data_UQPH;
2892  * ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
2893  *
2894  * @endcode
2895  *
2896  * We then pass them and the one-cell update function to the WorkStream to
2897  * be processed:
2898  *
2899  * @code
2900  * WorkStream::run(dof_handler.active_cell_iterators(),
2901  * *this,
2902  * &Solid::update_qph_incremental_one_cell,
2903  * &Solid::copy_local_to_global_UQPH,
2904  * scratch_data_UQPH,
2905  * per_task_data_UQPH);
2906  *
2907  * timer.leave_subsection();
2908  * }
2909  *
2910  *
2911  * @endcode
2912  *
2913  * Now we describe how we extract data from the solution vector and pass it
2914  * along to each QP storage object for processing.
2915  *
2916  * @code
2917  * template <int dim>
2918  * void Solid<dim>::update_qph_incremental_one_cell(
2919  * const typename DoFHandler<dim>::active_cell_iterator &cell,
2920  * ScratchData_UQPH & scratch,
2921  * PerTaskData_UQPH & /*data*/)
2922  * {
2923  * const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
2924  * quadrature_point_history.get_data(cell);
2925  * Assert(lqph.size() == n_q_points, ExcInternalError());
2926  *
2927  * Assert(scratch.solution_grads_u_total.size() == n_q_points,
2928  * ExcInternalError());
2929  * Assert(scratch.solution_values_p_total.size() == n_q_points,
2930  * ExcInternalError());
2931  * Assert(scratch.solution_values_J_total.size() == n_q_points,
2932  * ExcInternalError());
2933  *
2934  * scratch.reset();
2935  *
2936  * @endcode
2937  *
2938  * We first need to find the values and gradients at quadrature points
2939  * inside the current cell and then we update each local QP using the
2940  * displacement gradient and total pressure and dilatation solution
2941  * values:
2942  *
2943  * @code
2944  * scratch.fe_values.reinit(cell);
2945  * scratch.fe_values[u_fe].get_function_gradients(
2946  * scratch.solution_total, scratch.solution_grads_u_total);
2947  * scratch.fe_values[p_fe].get_function_values(
2948  * scratch.solution_total, scratch.solution_values_p_total);
2949  * scratch.fe_values[J_fe].get_function_values(
2950  * scratch.solution_total, scratch.solution_values_J_total);
2951  *
2952  * for (const unsigned int q_point :
2953  * scratch.fe_values.quadrature_point_indices())
2954  * lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point],
2955  * scratch.solution_values_p_total[q_point],
2956  * scratch.solution_values_J_total[q_point]);
2957  * }
2958  *
2959  *
2960  * @endcode
2961  *
2962  *
2963  * <a name="Solidsolve_nonlinear_timestep"></a>
2964  * <h4>Solid::solve_nonlinear_timestep</h4>
2965  *
2966 
2967  *
2968  * The next function is the driver method for the Newton-Raphson scheme. At
2969  * its top we create a new vector to store the current Newton update step,
2970  * reset the error storage objects and print solver header.
2971  *
2972  * @code
2973  * template <int dim>
2974  * void Solid<dim>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
2975  * {
2976  * std::cout << std::endl
2977  * << "Timestep " << time.get_timestep() << " @ " << time.current()
2978  * << "s" << std::endl;
2979  *
2980  * BlockVector<double> newton_update(dofs_per_block);
2981  *
2982  * error_residual.reset();
2983  * error_residual_0.reset();
2984  * error_residual_norm.reset();
2985  * error_update.reset();
2986  * error_update_0.reset();
2987  * error_update_norm.reset();
2988  *
2989  * print_conv_header();
2990  *
2991  * @endcode
2992  *
2993  * We now perform a number of Newton iterations to iteratively solve the
2994  * nonlinear problem. Since the problem is fully nonlinear and we are
2995  * using a full Newton method, the data stored in the tangent matrix and
2996  * right-hand side vector is not reusable and must be cleared at each
2997  * Newton step. We then initially build the right-hand side vector to
2998  * check for convergence (and store this value in the first iteration).
2999  * The unconstrained DOFs of the rhs vector hold the out-of-balance
3000  * forces. The building is done before assembling the system matrix as the
3001  * latter is an expensive operation and we can potentially avoid an extra
3002  * assembly process by not assembling the tangent matrix when convergence
3003  * is attained.
3004  *
3005  * @code
3006  * unsigned int newton_iteration = 0;
3007  * for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration)
3008  * {
3009  * std::cout << " " << std::setw(2) << newton_iteration << " "
3010  * << std::flush;
3011  *
3012  * tangent_matrix = 0.0;
3013  * system_rhs = 0.0;
3014  *
3015  * assemble_system_rhs();
3016  * get_error_residual(error_residual);
3017  *
3018  * if (newton_iteration == 0)
3019  * error_residual_0 = error_residual;
3020  *
3021  * @endcode
3022  *
3023  * We can now determine the normalized residual error and check for
3024  * solution convergence:
3025  *
3026  * @code
3027  * error_residual_norm = error_residual;
3028  * error_residual_norm.normalize(error_residual_0);
3029  *
3030  * if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u &&
3031  * error_residual_norm.u <= parameters.tol_f)
3032  * {
3033  * std::cout << " CONVERGED! " << std::endl;
3034  * print_conv_footer();
3035  *
3036  * break;
3037  * }
3038  *
3039  * @endcode
3040  *
3041  * If we have decided that we want to continue with the iteration, we
3042  * assemble the tangent, make and impose the Dirichlet constraints,
3043  * and do the solve of the linearised system:
3044  *
3045  * @code
3046  * assemble_system_tangent();
3047  * make_constraints(newton_iteration);
3048  * constraints.condense(tangent_matrix, system_rhs);
3049  *
3050  * const std::pair<unsigned int, double> lin_solver_output =
3051  * solve_linear_system(newton_update);
3052  *
3053  * get_error_update(newton_update, error_update);
3054  * if (newton_iteration == 0)
3055  * error_update_0 = error_update;
3056  *
3057  * @endcode
3058  *
3059  * We can now determine the normalized Newton update error, and
3060  * perform the actual update of the solution increment for the current
3061  * time step, update all quadrature point information pertaining to
3062  * this new displacement and stress state and continue iterating:
3063  *
3064  * @code
3065  * error_update_norm = error_update;
3066  * error_update_norm.normalize(error_update_0);
3067  *
3068  * solution_delta += newton_update;
3069  * update_qph_incremental(solution_delta);
3070  *
3071  * std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
3072  * << std::scientific << lin_solver_output.first << " "
3073  * << lin_solver_output.second << " "
3074  * << error_residual_norm.norm << " " << error_residual_norm.u
3075  * << " " << error_residual_norm.p << " "
3076  * << error_residual_norm.J << " " << error_update_norm.norm
3077  * << " " << error_update_norm.u << " " << error_update_norm.p
3078  * << " " << error_update_norm.J << " " << std::endl;
3079  * }
3080  *
3081  * @endcode
3082  *
3083  * At the end, if it turns out that we have in fact done more iterations
3084  * than the parameter file allowed, we raise an exception that can be
3085  * caught in the main() function. The call <code>AssertThrow(condition,
3086  * exc_object)</code> is in essence equivalent to <code>if (!cond) throw
3087  * exc_object;</code> but the former form fills certain fields in the
3088  * exception object that identify the location (filename and line number)
3089  * where the exception was raised to make it simpler to identify where the
3090  * problem happened.
3091  *
3092  * @code
3093  * AssertThrow(newton_iteration < parameters.max_iterations_NR,
3094  * ExcMessage("No convergence in nonlinear solver!"));
3095  * }
3096  *
3097  *
3098  * @endcode
3099  *
3100  *
3101  * <a name="Solidprint_conv_headerandSolidprint_conv_footer"></a>
3102  * <h4>Solid::print_conv_header and Solid::print_conv_footer</h4>
3103  *
3104 
3105  *
3106  * This program prints out data in a nice table that is updated
3107  * on a per-iteration basis. The next two functions set up the table
3108  * header and footer:
3109  *
3110  * @code
3111  * template <int dim>
3112  * void Solid<dim>::print_conv_header()
3113  * {
3114  * static const unsigned int l_width = 155;
3115  *
3116  * for (unsigned int i = 0; i < l_width; ++i)
3117  * std::cout << "_";
3118  * std::cout << std::endl;
3119  *
3120  * std::cout << " SOLVER STEP "
3121  * << " | LIN_IT LIN_RES RES_NORM "
3122  * << " RES_U RES_P RES_J NU_NORM "
3123  * << " NU_U NU_P NU_J " << std::endl;
3124  *
3125  * for (unsigned int i = 0; i < l_width; ++i)
3126  * std::cout << "_";
3127  * std::cout << std::endl;
3128  * }
3129  *
3130  *
3131  *
3132  * template <int dim>
3133  * void Solid<dim>::print_conv_footer()
3134  * {
3135  * static const unsigned int l_width = 155;
3136  *
3137  * for (unsigned int i = 0; i < l_width; ++i)
3138  * std::cout << "_";
3139  * std::cout << std::endl;
3140  *
3141  * const std::pair<double, double> error_dil = get_error_dilation();
3142  *
3143  * std::cout << "Relative errors:" << std::endl
3144  * << "Displacement:\t" << error_update.u / error_update_0.u
3145  * << std::endl
3146  * << "Force: \t\t" << error_residual.u / error_residual_0.u
3147  * << std::endl
3148  * << "Dilatation:\t" << error_dil.first << std::endl
3149  * << "v / V_0:\t" << error_dil.second * vol_reference << " / "
3150  * << vol_reference << " = " << error_dil.second << std::endl;
3151  * }
3152  *
3153  *
3154  * @endcode
3155  *
3156  *
3157  * <a name="Solidget_error_dilation"></a>
3158  * <h4>Solid::get_error_dilation</h4>
3159  *
3160 
3161  *
3162  * Calculate the volume of the domain in the spatial configuration
3163  *
3164  * @code
3165  * template <int dim>
3166  * double Solid<dim>::compute_vol_current() const
3167  * {
3168  * double vol_current = 0.0;
3169  *
3170  * FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
3171  *
3172  * for (const auto &cell : triangulation.active_cell_iterators())
3173  * {
3174  * fe_values.reinit(cell);
3175  *
3176  * @endcode
3177  *
3178  * In contrast to that which was previously called for,
3179  * in this instance the quadrature point data is specifically
3180  * non-modifiable since we will only be accessing data.
3181  * We ensure that the right get_data function is called by
3182  * marking this update function as constant.
3183  *
3184  * @code
3185  * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3186  * quadrature_point_history.get_data(cell);
3187  * Assert(lqph.size() == n_q_points, ExcInternalError());
3188  *
3189  * for (const unsigned int q_point : fe_values.quadrature_point_indices())
3190  * {
3191  * const double det_F_qp = lqph[q_point]->get_det_F();
3192  * const double JxW = fe_values.JxW(q_point);
3193  *
3194  * vol_current += det_F_qp * JxW;
3195  * }
3196  * }
3197  * Assert(vol_current > 0.0, ExcInternalError());
3198  * return vol_current;
3199  * }
3200  *
3201  * @endcode
3202  *
3203  * Calculate how well the dilatation @f$\widetilde{J}@f$ agrees with @f$J
3204  * \dealcoloneq \textrm{det}\ \mathbf{F}@f$ from the @f$L^2@f$ error @f$ \bigl[
3205  * \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}@f$.
3206  * We also return the ratio of the current volume of the
3207  * domain to the reference volume. This is of interest for incompressible
3208  * media where we want to check how well the isochoric constraint has been
3209  * enforced.
3210  *
3211  * @code
3212  * template <int dim>
3213  * std::pair<double, double> Solid<dim>::get_error_dilation() const
3214  * {
3215  * double dil_L2_error = 0.0;
3216  *
3217  * FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
3218  *
3219  * for (const auto &cell : triangulation.active_cell_iterators())
3220  * {
3221  * fe_values.reinit(cell);
3222  *
3223  * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3224  * quadrature_point_history.get_data(cell);
3225  * Assert(lqph.size() == n_q_points, ExcInternalError());
3226  *
3227  * for (const unsigned int q_point : fe_values.quadrature_point_indices())
3228  * {
3229  * const double det_F_qp = lqph[q_point]->get_det_F();
3230  * const double J_tilde_qp = lqph[q_point]->get_J_tilde();
3231  * const double the_error_qp_squared =
3232  * std::pow((det_F_qp - J_tilde_qp), 2);
3233  * const double JxW = fe_values.JxW(q_point);
3234  *
3235  * dil_L2_error += the_error_qp_squared * JxW;
3236  * }
3237  * }
3238  *
3239  * return std::make_pair(std::sqrt(dil_L2_error),
3240  * compute_vol_current() / vol_reference);
3241  * }
3242  *
3243  *
3244  * @endcode
3245  *
3246  *
3247  * <a name="Solidget_error_residual"></a>
3248  * <h4>Solid::get_error_residual</h4>
3249  *
3250 
3251  *
3252  * Determine the true residual error for the problem. That is, determine the
3253  * error in the residual for the unconstrained degrees of freedom. Note that
3254  * to do so, we need to ignore constrained DOFs by setting the residual in
3255  * these vector components to zero.
3256  *
3257  * @code
3258  * template <int dim>
3259  * void Solid<dim>::get_error_residual(Errors &error_residual)
3260  * {
3261  * BlockVector<double> error_res(dofs_per_block);
3262  *
3263  * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
3264  * if (!constraints.is_constrained(i))
3265  * error_res(i) = system_rhs(i);
3266  *
3267  * error_residual.norm = error_res.l2_norm();
3268  * error_residual.u = error_res.block(u_dof).l2_norm();
3269  * error_residual.p = error_res.block(p_dof).l2_norm();
3270  * error_residual.J = error_res.block(J_dof).l2_norm();
3271  * }
3272  *
3273  *
3274  * @endcode
3275  *
3276  *
3277  * <a name="Solidget_error_update"></a>
3278  * <h4>Solid::get_error_update</h4>
3279  *
3280 
3281  *
3282  * Determine the true Newton update error for the problem
3283  *
3284  * @code
3285  * template <int dim>
3286  * void Solid<dim>::get_error_update(const BlockVector<double> &newton_update,
3287  * Errors & error_update)
3288  * {
3289  * BlockVector<double> error_ud(dofs_per_block);
3290  * for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
3291  * if (!constraints.is_constrained(i))
3292  * error_ud(i) = newton_update(i);
3293  *
3294  * error_update.norm = error_ud.l2_norm();
3295  * error_update.u = error_ud.block(u_dof).l2_norm();
3296  * error_update.p = error_ud.block(p_dof).l2_norm();
3297  * error_update.J = error_ud.block(J_dof).l2_norm();
3298  * }
3299  *
3300  *
3301  *
3302  * @endcode
3303  *
3304  *
3305  * <a name="Solidget_total_solution"></a>
3306  * <h4>Solid::get_total_solution</h4>
3307  *
3308 
3309  *
3310  * This function provides the total solution, which is valid at any Newton
3311  * step. This is required as, to reduce computational error, the total
3312  * solution is only updated at the end of the timestep.
3313  *
3314  * @code
3315  * template <int dim>
3316  * BlockVector<double> Solid<dim>::get_total_solution(
3317  * const BlockVector<double> &solution_delta) const
3318  * {
3319  * BlockVector<double> solution_total(solution_n);
3320  * solution_total += solution_delta;
3321  * return solution_total;
3322  * }
3323  *
3324  *
3325  * @endcode
3326  *
3327  *
3328  * <a name="Solidassemble_system_tangent"></a>
3329  * <h4>Solid::assemble_system_tangent</h4>
3330  *
3331 
3332  *
3333  * Since we use TBB for assembly, we simply setup a copy of the
3334  * data structures required for the process and pass them, along
3335  * with the memory addresses of the assembly functions to the
3336  * WorkStream object for processing. Note that we must ensure that
3337  * the matrix is reset before any assembly operations can occur.
3338  *
3339  * @code
3340  * template <int dim>
3341  * void Solid<dim>::assemble_system_tangent()
3342  * {
3343  * timer.enter_subsection("Assemble tangent matrix");
3344  * std::cout << " ASM_K " << std::flush;
3345  *
3346  * tangent_matrix = 0.0;
3347  *
3348  * const UpdateFlags uf_cell(update_values | update_gradients |
3349  * update_JxW_values);
3350  *
3351  * PerTaskData_K per_task_data(dofs_per_cell);
3352  * ScratchData_K scratch_data(fe, qf_cell, uf_cell);
3353  *
3354  * @endcode
3355  *
3356  * The syntax used here to pass data to the WorkStream class
3357  * is discussed in @ref step_14 "step-14". We need to use this particular
3358  * call to WorkStream because assemble_system_tangent_one_cell
3359  * is a constant function and copy_local_to_global_K is
3360  * non-constant.
3361  *
3362  * @code
3363  * WorkStream::run(
3364  * dof_handler.active_cell_iterators(),
3365  * [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
3366  * ScratchData_K & scratch,
3367  * PerTaskData_K & data) {
3368  * this->assemble_system_tangent_one_cell(cell, scratch, data);
3369  * },
3370  * [this](const PerTaskData_K &data) { this->copy_local_to_global_K(data); },
3371  * scratch_data,
3372  * per_task_data);
3373  *
3374  * timer.leave_subsection();
3375  * }
3376  *
3377  * @endcode
3378  *
3379  * This function adds the local contribution to the system matrix.
3380  * Note that we choose not to use the constraint matrix to do the
3381  * job for us because the tangent matrix and residual processes have
3382  * been split up into two separate functions.
3383  *
3384  * @code
3385  * template <int dim>
3386  * void Solid<dim>::copy_local_to_global_K(const PerTaskData_K &data)
3387  * {
3388  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
3389  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
3390  * tangent_matrix.add(data.local_dof_indices[i],
3391  * data.local_dof_indices[j],
3392  * data.cell_matrix(i, j));
3393  * }
3394  *
3395  * @endcode
3396  *
3397  * Of course, we still have to define how we assemble the tangent matrix
3398  * contribution for a single cell. We first need to reset and initialize some
3399  * of the scratch data structures and retrieve some basic information
3400  * regarding the DOF numbering on this cell. We can precalculate the cell
3401  * shape function values and gradients. Note that the shape function gradients
3402  * are defined with regard to the current configuration. That is
3403  * @f$\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi}
3404  * \ \mathbf{F}^{-1}@f$.
3405  *
3406  * @code
3407  * template <int dim>
3408  * void Solid<dim>::assemble_system_tangent_one_cell(
3409  * const typename DoFHandler<dim>::active_cell_iterator &cell,
3410  * ScratchData_K & scratch,
3411  * PerTaskData_K & data) const
3412  * {
3413  * data.reset();
3414  * scratch.reset();
3415  * scratch.fe_values.reinit(cell);
3416  * cell->get_dof_indices(data.local_dof_indices);
3417  *
3418  * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3419  * quadrature_point_history.get_data(cell);
3420  * Assert(lqph.size() == n_q_points, ExcInternalError());
3421  *
3422  * for (const unsigned int q_point :
3423  * scratch.fe_values.quadrature_point_indices())
3424  * {
3425  * const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
3426  * for (const unsigned int k : scratch.fe_values.dof_indices())
3427  * {
3428  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
3429  *
3430  * if (k_group == u_dof)
3431  * {
3432  * scratch.grad_Nx[q_point][k] =
3433  * scratch.fe_values[u_fe].gradient(k, q_point) * F_inv;
3434  * scratch.symm_grad_Nx[q_point][k] =
3435  * symmetrize(scratch.grad_Nx[q_point][k]);
3436  * }
3437  * else if (k_group == p_dof)
3438  * scratch.Nx[q_point][k] =
3439  * scratch.fe_values[p_fe].value(k, q_point);
3440  * else if (k_group == J_dof)
3441  * scratch.Nx[q_point][k] =
3442  * scratch.fe_values[J_fe].value(k, q_point);
3443  * else
3444  * Assert(k_group <= J_dof, ExcInternalError());
3445  * }
3446  * }
3447  *
3448  * @endcode
3449  *
3450  * Now we build the local cell stiffness matrix. Since the global and
3451  * local system matrices are symmetric, we can exploit this property by
3452  * building only the lower half of the local matrix and copying the values
3453  * to the upper half. So we only assemble half of the
3454  * @f$\mathsf{\mathbf{k}}_{uu}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
3455  * \widetilde{p}} = \mathbf{0}@f$, @f$\mathsf{\mathbf{k}}_{\widetilde{J}
3456  * \widetilde{J}}@f$ blocks, while the whole
3457  * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$,
3458  * @f$\mathsf{\mathbf{k}}_{u \widetilde{J}} = \mathbf{0}@f$,
3459  * @f$\mathsf{\mathbf{k}}_{u \widetilde{p}}@f$ blocks are built.
3460  *
3461 
3462  *
3463  * In doing so, we first extract some configuration dependent variables
3464  * from our quadrature history objects for the current quadrature point.
3465  *
3466  * @code
3467  * for (const unsigned int q_point :
3468  * scratch.fe_values.quadrature_point_indices())
3469  * {
3470  * const Tensor<2, dim> tau = lqph[q_point]->get_tau();
3471  * const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc();
3472  * const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2();
3473  * const double det_F = lqph[q_point]->get_det_F();
3474  * const SymmetricTensor<2, dim> &I =
3475  * Physics::Elasticity::StandardTensors<dim>::I;
3476  *
3477  * @endcode
3478  *
3479  * Next we define some aliases to make the assembly process easier to
3480  * follow
3481  *
3482  * @code
3483  * const std::vector<double> & N = scratch.Nx[q_point];
3484  * const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
3485  * scratch.symm_grad_Nx[q_point];
3486  * const std::vector<Tensor<2, dim>> &grad_Nx = scratch.grad_Nx[q_point];
3487  * const double JxW = scratch.fe_values.JxW(q_point);
3488  *
3489  * for (const unsigned int i : scratch.fe_values.dof_indices())
3490  * {
3491  * const unsigned int component_i =
3492  * fe.system_to_component_index(i).first;
3493  * const unsigned int i_group = fe.system_to_base_index(i).first.first;
3494  *
3495  * for (const unsigned int j :
3496  * scratch.fe_values.dof_indices_ending_at(i))
3497  * {
3498  * const unsigned int component_j =
3499  * fe.system_to_component_index(j).first;
3500  * const unsigned int j_group =
3501  * fe.system_to_base_index(j).first.first;
3502  *
3503  * @endcode
3504  *
3505  * This is the @f$\mathsf{\mathbf{k}}_{uu}@f$
3506  * contribution. It comprises a material contribution, and a
3507  * geometrical stress contribution which is only added along
3508  * the local matrix diagonals:
3509  *
3510  * @code
3511  * if ((i_group == j_group) && (i_group == u_dof))
3512  * {
3513  * @endcode
3514  *
3515  * The material contribution:
3516  *
3517  * @code
3518  * data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc *
3519  * symm_grad_Nx[j] * JxW;
3520  *
3521  * @endcode
3522  *
3523  * The geometrical stress contribution:
3524  *
3525  * @code
3526  * if (component_i == component_j)
3527  * data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau *
3528  * grad_Nx[j][component_j] * JxW;
3529  * }
3530  * @endcode
3531  *
3532  * Next is the @f$\mathsf{\mathbf{k}}_{ \widetilde{p} u}@f$
3533  * contribution
3534  *
3535  * @code
3536  * else if ((i_group == p_dof) && (j_group == u_dof))
3537  * {
3538  * data.cell_matrix(i, j) += N[i] * det_F *
3539  * (symm_grad_Nx[j] * I) * JxW;
3540  * }
3541  * @endcode
3542  *
3543  * and lastly the @f$\mathsf{\mathbf{k}}_{ \widetilde{J}
3544  * \widetilde{p}}@f$ and @f$\mathsf{\mathbf{k}}_{ \widetilde{J}
3545  * \widetilde{J}}@f$ contributions:
3546  *
3547  * @code
3548  * else if ((i_group == J_dof) && (j_group == p_dof))
3549  * data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
3550  * else if ((i_group == j_group) && (i_group == J_dof))
3551  * data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
3552  * else
3553  * Assert((i_group <= J_dof) && (j_group <= J_dof),
3554  * ExcInternalError());
3555  * }
3556  * }
3557  * }
3558  *
3559  * @endcode
3560  *
3561  * Finally, we need to copy the lower half of the local matrix into the
3562  * upper half:
3563  *
3564  * @code
3565  * for (const unsigned int i : scratch.fe_values.dof_indices())
3566  * for (const unsigned int j :
3567  * scratch.fe_values.dof_indices_starting_at(i + 1))
3568  * data.cell_matrix(i, j) = data.cell_matrix(j, i);
3569  * }
3570  *
3571  * @endcode
3572  *
3573  *
3574  * <a name="Solidassemble_system_rhs"></a>
3575  * <h4>Solid::assemble_system_rhs</h4>
3576  * The assembly of the right-hand side process is similar to the
3577  * tangent matrix, so we will not describe it in too much detail.
3578  * Note that since we are describing a problem with Neumann BCs,
3579  * we will need the face normals and so must specify this in the
3580  * update flags.
3581  *
3582  * @code
3583  * template <int dim>
3584  * void Solid<dim>::assemble_system_rhs()
3585  * {
3586  * timer.enter_subsection("Assemble system right-hand side");
3587  * std::cout << " ASM_R " << std::flush;
3588  *
3589  * system_rhs = 0.0;
3590  *
3591  * const UpdateFlags uf_cell(update_values | update_gradients |
3592  * update_JxW_values);
3593  * const UpdateFlags uf_face(update_values | update_normal_vectors |
3594  * update_JxW_values);
3595  *
3596  * PerTaskData_RHS per_task_data(dofs_per_cell);
3597  * ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
3598  *
3599  * WorkStream::run(
3600  * dof_handler.active_cell_iterators(),
3601  * [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
3602  * ScratchData_RHS & scratch,
3603  * PerTaskData_RHS & data) {
3604  * this->assemble_system_rhs_one_cell(cell, scratch, data);
3605  * },
3606  * [this](const PerTaskData_RHS &data) {
3607  * this->copy_local_to_global_rhs(data);
3608  * },
3609  * scratch_data,
3610  * per_task_data);
3611  *
3612  * timer.leave_subsection();
3613  * }
3614  *
3615  *
3616  *
3617  * template <int dim>
3618  * void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS &data)
3619  * {
3620  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
3621  * system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
3622  * }
3623  *
3624  *
3625  *
3626  * template <int dim>
3627  * void Solid<dim>::assemble_system_rhs_one_cell(
3628  * const typename DoFHandler<dim>::active_cell_iterator &cell,
3629  * ScratchData_RHS & scratch,
3630  * PerTaskData_RHS & data) const
3631  * {
3632  * data.reset();
3633  * scratch.reset();
3634  * scratch.fe_values.reinit(cell);
3635  * cell->get_dof_indices(data.local_dof_indices);
3636  *
3637  * const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
3638  * quadrature_point_history.get_data(cell);
3639  * Assert(lqph.size() == n_q_points, ExcInternalError());
3640  *
3641  * for (const unsigned int q_point :
3642  * scratch.fe_values.quadrature_point_indices())
3643  * {
3644  * const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
3645  *
3646  * for (const unsigned int k : scratch.fe_values.dof_indices())
3647  * {
3648  * const unsigned int k_group = fe.system_to_base_index(k).first.first;
3649  *
3650  * if (k_group == u_dof)
3651  * scratch.symm_grad_Nx[q_point][k] = symmetrize(
3652  * scratch.fe_values[u_fe].gradient(k, q_point) * F_inv);
3653  * else if (k_group == p_dof)
3654  * scratch.Nx[q_point][k] =
3655  * scratch.fe_values[p_fe].value(k, q_point);
3656  * else if (k_group == J_dof)
3657  * scratch.Nx[q_point][k] =
3658  * scratch.fe_values[J_fe].value(k, q_point);
3659  * else
3660  * Assert(k_group <= J_dof, ExcInternalError());
3661  * }
3662  * }
3663  *
3664  * for (const unsigned int q_point :
3665  * scratch.fe_values.quadrature_point_indices())
3666  * {
3667  * const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau();
3668  * const double det_F = lqph[q_point]->get_det_F();
3669  * const double J_tilde = lqph[q_point]->get_J_tilde();
3670  * const double p_tilde = lqph[q_point]->get_p_tilde();
3671  * const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ();
3672  *
3673  * const std::vector<double> & N = scratch.Nx[q_point];
3674  * const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
3675  * scratch.symm_grad_Nx[q_point];
3676  * const double JxW = scratch.fe_values.JxW(q_point);
3677  *
3678  * @endcode
3679  *
3680  * We first compute the contributions
3681  * from the internal forces. Note, by
3682  * definition of the rhs as the negative
3683  * of the residual, these contributions
3684  * are subtracted.
3685  *
3686  * @code
3687  * for (const unsigned int i : scratch.fe_values.dof_indices())
3688  * {
3689  * const unsigned int i_group = fe.system_to_base_index(i).first.first;
3690  *
3691  * if (i_group == u_dof)
3692  * data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
3693  * else if (i_group == p_dof)
3694  * data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
3695  * else if (i_group == J_dof)
3696  * data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW;
3697  * else
3698  * Assert(i_group <= J_dof, ExcInternalError());
3699  * }
3700  * }
3701  *
3702  * @endcode
3703  *
3704  * Next we assemble the Neumann contribution. We first check to see it the
3705  * cell face exists on a boundary on which a traction is applied and add
3706  * the contribution if this is the case.
3707  *
3708  * @code
3709  * for (const auto &face : cell->face_iterators())
3710  * if (face->at_boundary() && face->boundary_id() == 6)
3711  * {
3712  * scratch.fe_face_values.reinit(cell, face);
3713  *
3714  * for (const unsigned int f_q_point :
3715  * scratch.fe_face_values.quadrature_point_indices())
3716  * {
3717  * const Tensor<1, dim> &N =
3718  * scratch.fe_face_values.normal_vector(f_q_point);
3719  *
3720  * @endcode
3721  *
3722  * Using the face normal at this quadrature point we specify the
3723  * traction in reference configuration. For this problem, a
3724  * defined pressure is applied in the reference configuration.
3725  * The direction of the applied traction is assumed not to
3726  * evolve with the deformation of the domain. The traction is
3727  * defined using the first Piola-Kirchhoff stress is simply
3728  * @f$\mathbf{t} = \mathbf{P}\mathbf{N} = [p_0 \mathbf{I}]
3729  * \mathbf{N} = p_0 \mathbf{N}@f$ We use the time variable to
3730  * linearly ramp up the pressure load.
3731  *
3732 
3733  *
3734  * Note that the contributions to the right hand side vector we
3735  * compute here only exist in the displacement components of the
3736  * vector.
3737  *
3738  * @code
3739  * static const double p0 =
3740  * -4.0 / (parameters.scale * parameters.scale);
3741  * const double time_ramp = (time.current() / time.end());
3742  * const double pressure = p0 * parameters.p_p0 * time_ramp;
3743  * const Tensor<1, dim> traction = pressure * N;
3744  *
3745  * for (const unsigned int i : scratch.fe_values.dof_indices())
3746  * {
3747  * const unsigned int i_group =
3748  * fe.system_to_base_index(i).first.first;
3749  *
3750  * if (i_group == u_dof)
3751  * {
3752  * const unsigned int component_i =
3753  * fe.system_to_component_index(i).first;
3754  * const double Ni =
3755  * scratch.fe_face_values.shape_value(i, f_q_point);
3756  * const double JxW = scratch.fe_face_values.JxW(f_q_point);
3757  *
3758  * data.cell_rhs(i) += (Ni * traction[component_i]) * JxW;
3759  * }
3760  * }
3761  * }
3762  * }
3763  * }
3764  *
3765  * @endcode
3766  *
3767  *
3768  * <a name="Solidmake_constraints"></a>
3769  * <h4>Solid::make_constraints</h4>
3770  * The constraints for this problem are simple to describe.
3771  * However, since we are dealing with an iterative Newton method,
3772  * it should be noted that any displacement constraints should only
3773  * be specified at the zeroth iteration and subsequently no
3774  * additional contributions are to be made since the constraints
3775  * are already exactly satisfied.
3776  *
3777  * @code
3778  * template <int dim>
3779  * void Solid<dim>::make_constraints(const int &it_nr)
3780  * {
3781  * std::cout << " CST " << std::flush;
3782  *
3783  * @endcode
3784  *
3785  * Since the constraints are different at different Newton iterations, we
3786  * need to clear the constraints matrix and completely rebuild
3787  * it. However, after the first iteration, the constraints remain the same
3788  * and we can simply skip the rebuilding step if we do not clear it.
3789  *
3790  * @code
3791  * if (it_nr > 1)
3792  * return;
3793  * constraints.clear();
3794  * const bool apply_dirichlet_bc = (it_nr == 0);
3795  *
3796  * @endcode
3797  *
3798  * In this particular example, the boundary values will be calculated for
3799  * the two first iterations of Newton's algorithm. In general, one would
3800  * build non-homogeneous constraints in the zeroth iteration (that is, when
3801  * `apply_dirichlet_bc == true`) and build only the corresponding
3802  * homogeneous constraints in the following step. While the current
3803  * example has only homogeneous constraints, previous experiences have
3804  * shown that a common error is forgetting to add the extra condition when
3805  * refactoring the code to specific uses. This could lead errors that are
3806  * hard to debug. In this spirit, we choose to make the code more verbose
3807  * in terms of what operations are performed at each Newton step.
3808  *
3809 
3810  *
3811  * The boundary conditions for the indentation problem are as follows: On
3812  * the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry condition to
3813  * allow only planar movement while the +x and +z faces (IDs 1,5) are
3814  * traction free. In this contrived problem, part of the +y face (ID 3) is
3815  * set to have no motion in the x- and z-component. Finally, as described
3816  * earlier, the other part of the +y face has an the applied pressure but
3817  * is also constrained in the x- and z-directions.
3818  *
3819 
3820  *
3821  * In the following, we will have to tell the function interpolation
3822  * boundary values which components of the solution vector should be
3823  * constrained (i.e., whether it's the x-, y-, z-displacements or
3824  * combinations thereof). This is done using ComponentMask objects (see
3825  * @ref GlossComponentMask) which we can get from the finite element if we
3826  * provide it with an extractor object for the component we wish to
3827  * select. To this end we first set up such extractor objects and later
3828  * use it when generating the relevant component masks:
3829  *
3830  * @code
3831  * const FEValuesExtractors::Scalar x_displacement(0);
3832  * const FEValuesExtractors::Scalar y_displacement(1);
3833  *
3834  * {
3835  * const int boundary_id = 0;
3836  *
3837  * if (apply_dirichlet_bc == true)
3839  * dof_handler,
3840  * boundary_id,
3842  * constraints,
3843  * fe.component_mask(x_displacement));
3844  * else
3846  * dof_handler,
3847  * boundary_id,
3849  * constraints,
3850  * fe.component_mask(x_displacement));
3851  * }
3852  * {
3853  * const int boundary_id = 2;
3854  *
3855  * if (apply_dirichlet_bc == true)
3857  * dof_handler,
3858  * boundary_id,
3860  * constraints,
3861  * fe.component_mask(y_displacement));
3862  * else
3864  * dof_handler,
3865  * boundary_id,
3867  * constraints,
3868  * fe.component_mask(y_displacement));
3869  * }
3870  *
3871  * if (dim == 3)
3872  * {
3873  * const FEValuesExtractors::Scalar z_displacement(2);
3874  *
3875  * {
3876  * const int boundary_id = 3;
3877  *
3878  * if (apply_dirichlet_bc == true)
3880  * dof_handler,
3881  * boundary_id,
3883  * constraints,
3884  * (fe.component_mask(x_displacement) |
3885  * fe.component_mask(z_displacement)));
3886  * else
3888  * dof_handler,
3889  * boundary_id,
3891  * constraints,
3892  * (fe.component_mask(x_displacement) |
3893  * fe.component_mask(z_displacement)));
3894  * }
3895  * {
3896  * const int boundary_id = 4;
3897  *
3898  * if (apply_dirichlet_bc == true)
3900  * dof_handler,
3901  * boundary_id,
3903  * constraints,
3904  * fe.component_mask(z_displacement));
3905  * else
3907  * dof_handler,
3908  * boundary_id,
3910  * constraints,
3911  * fe.component_mask(z_displacement));
3912  * }
3913  *
3914  * {
3915  * const int boundary_id = 6;
3916  *
3917  * if (apply_dirichlet_bc == true)
3919  * dof_handler,
3920  * boundary_id,
3922  * constraints,
3923  * (fe.component_mask(x_displacement) |
3924  * fe.component_mask(z_displacement)));
3925  * else
3927  * dof_handler,
3928  * boundary_id,
3930  * constraints,
3931  * (fe.component_mask(x_displacement) |
3932  * fe.component_mask(z_displacement)));
3933  * }
3934  * }
3935  * else
3936  * {
3937  * {
3938  * const int boundary_id = 3;
3939  *
3940  * if (apply_dirichlet_bc == true)
3942  * dof_handler,
3943  * boundary_id,
3945  * constraints,
3946  * (fe.component_mask(x_displacement)));
3947  * else
3949  * dof_handler,
3950  * boundary_id,
3952  * constraints,
3953  * (fe.component_mask(x_displacement)));
3954  * }
3955  * {
3956  * const int boundary_id = 6;
3957  *
3958  * if (apply_dirichlet_bc == true)
3960  * dof_handler,
3961  * boundary_id,
3963  * constraints,
3964  * (fe.component_mask(x_displacement)));
3965  * else
3967  * dof_handler,
3968  * boundary_id,
3970  * constraints,
3971  * (fe.component_mask(x_displacement)));
3972  * }
3973  * }
3974  *
3975  * constraints.close();
3976  * }
3977  *
3978  * @endcode
3979  *
3980  *
3981  * <a name="Solidassemble_sc"></a>
3982  * <h4>Solid::assemble_sc</h4>
3983  * Solving the entire block system is a bit problematic as there are no
3984  * contributions to the @f$\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}@f$
3985  * block, rendering it noninvertible (when using an iterative solver).
3986  * Since the pressure and dilatation variables DOFs are discontinuous, we can
3987  * condense them out to form a smaller displacement-only system which
3988  * we will then solve and subsequently post-process to retrieve the
3989  * pressure and dilatation solutions.
3990  *
3991 
3992  *
3993  * The static condensation process could be performed at a global level but we
3994  * need the inverse of one of the blocks. However, since the pressure and
3995  * dilatation variables are discontinuous, the static condensation (SC)
3996  * operation can also be done on a per-cell basis and we can produce the
3997  * inverse of the block-diagonal
3998  * @f$\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}@f$ block by inverting the
3999  * local blocks. We can again use TBB to do this since each operation will be
4000  * independent of one another.
4001  *
4002 
4003  *
4004  * Using the TBB via the WorkStream class, we assemble the contributions to
4005  * form
4006  * @f$
4007  * \mathsf{\mathbf{K}}_{\textrm{con}}
4008  * = \bigl[ \mathsf{\mathbf{K}}_{uu} +
4009  * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
4010  * @f$
4011  * from each element's contributions. These contributions are then added to
4012  * the global stiffness matrix. Given this description, the following two
4013  * functions should be clear:
4014  *
4015  * @code
4016  * template <int dim>
4017  * void Solid<dim>::assemble_sc()
4018  * {
4019  * timer.enter_subsection("Perform static condensation");
4020  * std::cout << " ASM_SC " << std::flush;
4021  *
4022  * PerTaskData_SC per_task_data(dofs_per_cell,
4023  * element_indices_u.size(),
4024  * element_indices_p.size(),
4025  * element_indices_J.size());
4026  * ScratchData_SC scratch_data;
4027  *
4028  * WorkStream::run(dof_handler.active_cell_iterators(),
4029  * *this,
4030  * &Solid::assemble_sc_one_cell,
4031  * &Solid::copy_local_to_global_sc,
4032  * scratch_data,
4033  * per_task_data);
4034  *
4035  * timer.leave_subsection();
4036  * }
4037  *
4038  *
4039  * template <int dim>
4040  * void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC &data)
4041  * {
4042  * for (unsigned int i = 0; i < dofs_per_cell; ++i)
4043  * for (unsigned int j = 0; j < dofs_per_cell; ++j)
4044  * tangent_matrix.add(data.local_dof_indices[i],
4045  * data.local_dof_indices[j],
4046  * data.cell_matrix(i, j));
4047  * }
4048  *
4049  *
4050  * @endcode
4051  *
4052  * Now we describe the static condensation process. As per usual, we must
4053  * first find out which global numbers the degrees of freedom on this cell
4054  * have and reset some data structures:
4055  *
4056  * @code
4057  * template <int dim>
4058  * void Solid<dim>::assemble_sc_one_cell(
4059  * const typename DoFHandler<dim>::active_cell_iterator &cell,
4060  * ScratchData_SC & scratch,
4061  * PerTaskData_SC & data)
4062  * {
4063  * data.reset();
4064  * scratch.reset();
4065  * cell->get_dof_indices(data.local_dof_indices);
4066  *
4067  * @endcode
4068  *
4069  * We now extract the contribution of the dofs associated with the current
4070  * cell to the global stiffness matrix. The discontinuous nature of the
4071  * @f$\widetilde{p}@f$ and @f$\widetilde{J}@f$ interpolations mean that their is
4072  * no coupling of the local contributions at the global level. This is not
4073  * the case with the @f$\mathbf{u}@f$ dof. In other words,
4074  * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$,
4075  * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{p}}@f$ and
4076  * @f$\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}@f$, when extracted
4077  * from the global stiffness matrix are the element contributions. This
4078  * is not the case for @f$\mathsf{\mathbf{k}}_{uu}@f$.
4079  *
4080 
4081  *
4082  * Note: A lower-case symbol is used to denote element stiffness matrices.
4083  *
4084 
4085  *
4086  * Currently the matrix corresponding to
4087  * the dof associated with the current element
4088  * (denoted somewhat loosely as @f$\mathsf{\mathbf{k}}@f$)
4089  * is of the form:
4090  * @f{align*}
4091  * \begin{bmatrix}
4092  * \mathsf{\mathbf{k}}_{uu} & \mathsf{\mathbf{k}}_{u\widetilde{p}}
4093  * & \mathbf{0}
4094  * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
4095  * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}
4096  * \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
4097  * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
4098  * @f}
4099  *
4100 
4101  *
4102  * We now need to modify it such that it appear as
4103  * @f{align*}
4104  * \begin{bmatrix}
4105  * \mathsf{\mathbf{k}}_{\textrm{con}} &
4106  * \mathsf{\mathbf{k}}_{u\widetilde{p}} & \mathbf{0}
4107  * \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
4108  * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
4109  * \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
4110  * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
4111  * @f}
4112  * with @f$\mathsf{\mathbf{k}}_{\textrm{con}} = \bigl[
4113  * \mathsf{\mathbf{k}}_{uu} +\overline{\overline{\mathsf{\mathbf{k}}}}~
4114  * \bigr]@f$ where @f$ \overline{\overline{\mathsf{\mathbf{k}}}}
4115  * \dealcoloneq \mathsf{\mathbf{k}}_{u\widetilde{p}}
4116  * \overline{\mathsf{\mathbf{k}}} \mathsf{\mathbf{k}}_{\widetilde{p}u}
4117  * @f$
4118  * and
4119  * @f$
4120  * \overline{\mathsf{\mathbf{k}}} =
4121  * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}}^{-1}
4122  * \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}}
4123  * \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
4124  * @f$.
4125  *
4126 
4127  *
4128  * At this point, we need to take note of
4129  * the fact that global data already exists
4130  * in the @f$\mathsf{\mathbf{K}}_{uu}@f$,
4131  * @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
4132  * and
4133  * @f$\mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{p}}@f$
4134  * sub-blocks. So if we are to modify them, we must account for the data
4135  * that is already there (i.e. simply add to it or remove it if
4136  * necessary). Since the copy_local_to_global operation is a "+="
4137  * operation, we need to take this into account
4138  *
4139 
4140  *
4141  * For the @f$\mathsf{\mathbf{K}}_{uu}@f$ block in particular, this means that
4142  * contributions have been added from the surrounding cells, so we need to
4143  * be careful when we manipulate this block. We can't just erase the
4144  * sub-blocks.
4145  *
4146 
4147  *
4148  * This is the strategy we will employ to get the sub-blocks we want:
4149  *
4150 
4151  *
4152  * - @f$ {\mathsf{\mathbf{k}}}_{\textrm{store}}@f$:
4153  * Since we don't have access to @f$\mathsf{\mathbf{k}}_{uu}@f$,
4154  * but we know its contribution is added to
4155  * the global @f$\mathsf{\mathbf{K}}_{uu}@f$ matrix, we just want
4156  * to add the element wise
4157  * static-condensation @f$\overline{\overline{\mathsf{\mathbf{k}}}}@f$.
4158  *
4159 
4160  *
4161  * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$:
4162  * Similarly, @f$\mathsf{\mathbf{k}}_{\widetilde{p}
4163  * \widetilde{J}}@f$ exists in
4164  * the subblock. Since the copy
4165  * operation is a += operation, we
4166  * need to subtract the existing
4167  * @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
4168  * submatrix in addition to
4169  * "adding" that which we wish to
4170  * replace it with.
4171  *
4172 
4173  *
4174  * - @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}@f$:
4175  * Since the global matrix
4176  * is symmetric, this block is the
4177  * same as the one above and we
4178  * can simply use
4179  * @f$\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
4180  * as a substitute for this one.
4181  *
4182 
4183  *
4184  * We first extract element data from the
4185  * system matrix. So first we get the
4186  * entire subblock for the cell, then
4187  * extract @f$\mathsf{\mathbf{k}}@f$
4188  * for the dofs associated with
4189  * the current element
4190  *
4191  * @code
4192  * data.k_orig.extract_submatrix_from(tangent_matrix,
4193  * data.local_dof_indices,
4194  * data.local_dof_indices);
4195  * @endcode
4196  *
4197  * and next the local matrices for
4198  * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} u}@f$
4199  * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
4200  * and
4201  * @f$\mathsf{\mathbf{k}}_{ \widetilde{J} \widetilde{J}}@f$:
4202  *
4203  * @code
4204  * data.k_pu.extract_submatrix_from(data.k_orig,
4205  * element_indices_p,
4206  * element_indices_u);
4207  * data.k_pJ.extract_submatrix_from(data.k_orig,
4208  * element_indices_p,
4209  * element_indices_J);
4210  * data.k_JJ.extract_submatrix_from(data.k_orig,
4211  * element_indices_J,
4212  * element_indices_J);
4213  *
4214  * @endcode
4215  *
4216  * To get the inverse of @f$\mathsf{\mathbf{k}}_{\widetilde{p}
4217  * \widetilde{J}}@f$, we invert it directly. This operation is relatively
4218  * inexpensive since @f$\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}@f$
4219  * since block-diagonal.
4220  *
4221  * @code
4222  * data.k_pJ_inv.invert(data.k_pJ);
4223  *
4224  * @endcode
4225  *
4226  * Now we can make condensation terms to
4227  * add to the @f$\mathsf{\mathbf{k}}_{uu}@f$
4228  * block and put them in
4229  * the cell local matrix
4230  * @f$
4231  * \mathsf{\mathbf{A}}
4232  * =
4233  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4234  * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4235  * @f$:
4236  *
4237  * @code
4238  * data.k_pJ_inv.mmult(data.A, data.k_pu);
4239  * @endcode
4240  *
4241  * @f$
4242  * \mathsf{\mathbf{B}}
4243  * =
4244  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4245  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4246  * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4247  * @f$
4248  *
4249  * @code
4250  * data.k_JJ.mmult(data.B, data.A);
4251  * @endcode
4252  *
4253  * @f$
4254  * \mathsf{\mathbf{C}}
4255  * =
4256  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
4257  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4258  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4259  * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4260  * @f$
4261  *
4262  * @code
4263  * data.k_pJ_inv.Tmmult(data.C, data.B);
4264  * @endcode
4265  *
4266  * @f$
4267  * \overline{\overline{\mathsf{\mathbf{k}}}}
4268  * =
4269  * \mathsf{\mathbf{k}}_{u \widetilde{p}}
4270  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
4271  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
4272  * \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
4273  * \mathsf{\mathbf{k}}_{\widetilde{p} u}
4274  * @f$
4275  *
4276  * @code
4277  * data.k_pu.Tmmult(data.k_bbar, data.C);
4278  * data.k_bbar.scatter_matrix_to(element_indices_u,
4279  * element_indices_u,
4280  * data.cell_matrix);
4281  *
4282  * @endcode
4283  *
4284  * Next we place
4285  * @f$\mathsf{\mathbf{k}}^{-1}_{ \widetilde{p} \widetilde{J}}@f$
4286  * in the
4287  * @f$\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}@f$
4288  * block for post-processing. Note again
4289  * that we need to remove the
4290  * contribution that already exists there.
4291  *
4292  * @code
4293  * data.k_pJ_inv.add(-1.0, data.k_pJ);
4294  * data.k_pJ_inv.scatter_matrix_to(element_indices_p,
4295  * element_indices_J,
4296  * data.cell_matrix);
4297  * }
4298  *
4299  * @endcode
4300  *
4301  *
4302  * <a name="Solidsolve_linear_system"></a>
4303  * <h4>Solid::solve_linear_system</h4>
4304  * We now have all of the necessary components to use one of two possible
4305  * methods to solve the linearised system. The first is to perform static
4306  * condensation on an element level, which requires some alterations
4307  * to the tangent matrix and RHS vector. Alternatively, the full block
4308  * system can be solved by performing condensation on a global level.
4309  * Below we implement both approaches.
4310  *
4311  * @code
4312  * template <int dim>
4313  * std::pair<unsigned int, double>
4314  * Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
4315  * {
4316  * unsigned int lin_it = 0;
4317  * double lin_res = 0.0;
4318  *
4319  * if (parameters.use_static_condensation == true)
4320  * {
4321  * @endcode
4322  *
4323  * Firstly, here is the approach using the (permanent) augmentation of
4324  * the tangent matrix. For the following, recall that
4325  * @f{align*}
4326  * \mathsf{\mathbf{K}}_{\textrm{store}}
4327  * \dealcoloneq
4328  * \begin{bmatrix}
4329  * \mathsf{\mathbf{K}}_{\textrm{con}} &
4330  * \mathsf{\mathbf{K}}_{u\widetilde{p}} & \mathbf{0}
4331  * \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
4332  * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
4333  * \\ \mathbf{0} &
4334  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
4335  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, .
4336  * @f}
4337  * and
4338  * @f{align*}
4339  * d \widetilde{\mathsf{\mathbf{p}}}
4340  * & =
4341  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4342  * \bigl[
4343  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4344  * -
4345  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4346  * d \widetilde{\mathsf{\mathbf{J}}} \bigr]
4347  * \\ d \widetilde{\mathsf{\mathbf{J}}}
4348  * & =
4349  * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
4350  * \bigl[
4351  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4352  * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
4353  * \mathsf{\mathbf{u}} \bigr]
4354  * \\ \Rightarrow d \widetilde{\mathsf{\mathbf{p}}}
4355  * &= \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4356  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4357  * -
4358  * \underbrace{\bigl[\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4359  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4360  * \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}\bigr]}_{\overline{\mathsf{\mathbf{K}}}}\bigl[
4361  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4362  * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
4363  * \mathsf{\mathbf{u}} \bigr]
4364  * @f}
4365  * and thus
4366  * @f[
4367  * \underbrace{\bigl[ \mathsf{\mathbf{K}}_{uu} +
4368  * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
4369  * }_{\mathsf{\mathbf{K}}_{\textrm{con}}} d
4370  * \mathsf{\mathbf{u}}
4371  * =
4372  * \underbrace{
4373  * \Bigl[
4374  * \mathsf{\mathbf{F}}_{u}
4375  * - \mathsf{\mathbf{K}}_{u\widetilde{p}} \bigl[
4376  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4377  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4378  * -
4379  * \overline{\mathsf{\mathbf{K}}}\mathsf{\mathbf{F}}_{\widetilde{p}}
4380  * \bigr]
4381  * \Bigr]}_{\mathsf{\mathbf{F}}_{\textrm{con}}}
4382  * @f]
4383  * where
4384  * @f[
4385  * \overline{\overline{\mathsf{\mathbf{K}}}} \dealcoloneq
4386  * \mathsf{\mathbf{K}}_{u\widetilde{p}}
4387  * \overline{\mathsf{\mathbf{K}}}
4388  * \mathsf{\mathbf{K}}_{\widetilde{p}u} \, .
4389  * @f]
4390  *
4391 
4392  *
4393  * At the top, we allocate two temporary vectors to help with the
4394  * static condensation, and variables to store the number of
4395  * linear solver iterations and the (hopefully converged) residual.
4396  *
4397  * @code
4398  * BlockVector<double> A(dofs_per_block);
4399  * BlockVector<double> B(dofs_per_block);
4400  *
4401  *
4402  * @endcode
4403  *
4404  * In the first step of this function, we solve for the incremental
4405  * displacement @f$d\mathbf{u}@f$. To this end, we perform static
4406  * condensation to make
4407  * @f$\mathsf{\mathbf{K}}_{\textrm{con}}
4408  * = \bigl[ \mathsf{\mathbf{K}}_{uu} +
4409  * \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]@f$
4410  * and put
4411  * @f$\mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}@f$
4412  * in the original @f$\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}@f$
4413  * block. That is, we make @f$\mathsf{\mathbf{K}}_{\textrm{store}}@f$.
4414  *
4415  * @code
4416  * {
4417  * assemble_sc();
4418  *
4419  * @endcode
4420  *
4421  * @f$
4422  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4423  * =
4424  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4425  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4426  * @f$
4427  *
4428  * @code
4429  * tangent_matrix.block(p_dof, J_dof)
4430  * .vmult(A.block(J_dof), system_rhs.block(p_dof));
4431  * @endcode
4432  *
4433  * @f$
4434  * \mathsf{\mathbf{B}}_{\widetilde{J}}
4435  * =
4436  * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4437  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4438  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4439  * @f$
4440  *
4441  * @code
4442  * tangent_matrix.block(J_dof, J_dof)
4443  * .vmult(B.block(J_dof), A.block(J_dof));
4444  * @endcode
4445  *
4446  * @f$
4447  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4448  * =
4449  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4450  * -
4451  * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4452  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4453  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4454  * @f$
4455  *
4456  * @code
4457  * A.block(J_dof) = system_rhs.block(J_dof);
4458  * A.block(J_dof) -= B.block(J_dof);
4459  * @endcode
4460  *
4461  * @f$
4462  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4463  * =
4464  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4465  * [
4466  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4467  * -
4468  * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4469  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4470  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4471  * ]
4472  * @f$
4473  *
4474  * @code
4475  * tangent_matrix.block(p_dof, J_dof)
4476  * .Tvmult(A.block(p_dof), A.block(J_dof));
4477  * @endcode
4478  *
4479  * @f$
4480  * \mathsf{\mathbf{A}}_{u}
4481  * =
4482  * \mathsf{\mathbf{K}}_{u \widetilde{p}}
4483  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4484  * [
4485  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4486  * -
4487  * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4488  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4489  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4490  * ]
4491  * @f$
4492  *
4493  * @code
4494  * tangent_matrix.block(u_dof, p_dof)
4495  * .vmult(A.block(u_dof), A.block(p_dof));
4496  * @endcode
4497  *
4498  * @f$
4499  * \mathsf{\mathbf{F}}_{\text{con}}
4500  * =
4501  * \mathsf{\mathbf{F}}_{u}
4502  * -
4503  * \mathsf{\mathbf{K}}_{u \widetilde{p}}
4504  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
4505  * [
4506  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4507  * -
4508  * \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
4509  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
4510  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4511  * ]
4512  * @f$
4513  *
4514  * @code
4515  * system_rhs.block(u_dof) -= A.block(u_dof);
4516  *
4517  * timer.enter_subsection("Linear solver");
4518  * std::cout << " SLV " << std::flush;
4519  * if (parameters.type_lin == "CG")
4520  * {
4521  * const auto solver_its = static_cast<unsigned int>(
4522  * tangent_matrix.block(u_dof, u_dof).m() *
4523  * parameters.max_iterations_lin);
4524  * const double tol_sol =
4525  * parameters.tol_lin * system_rhs.block(u_dof).l2_norm();
4526  *
4527  * SolverControl solver_control(solver_its, tol_sol);
4528  *
4529  * GrowingVectorMemory<Vector<double>> GVM;
4530  * SolverCG<Vector<double>> solver_CG(solver_control, GVM);
4531  *
4532  * @endcode
4533  *
4534  * We've chosen by default a SSOR preconditioner as it appears to
4535  * provide the fastest solver convergence characteristics for this
4536  * problem on a single-thread machine. However, this might not be
4537  * true for different problem sizes.
4538  *
4539  * @code
4541  * preconditioner(parameters.preconditioner_type,
4542  * parameters.preconditioner_relaxation);
4543  * preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
4544  *
4545  * solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
4546  * newton_update.block(u_dof),
4547  * system_rhs.block(u_dof),
4548  * preconditioner);
4549  *
4550  * lin_it = solver_control.last_step();
4551  * lin_res = solver_control.last_value();
4552  * }
4553  * else if (parameters.type_lin == "Direct")
4554  * {
4555  * @endcode
4556  *
4557  * Otherwise if the problem is small
4558  * enough, a direct solver can be
4559  * utilised.
4560  *
4561  * @code
4562  * SparseDirectUMFPACK A_direct;
4563  * A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
4564  * A_direct.vmult(newton_update.block(u_dof),
4565  * system_rhs.block(u_dof));
4566  *
4567  * lin_it = 1;
4568  * lin_res = 0.0;
4569  * }
4570  * else
4571  * Assert(false, ExcMessage("Linear solver type not implemented"));
4572  *
4573  * timer.leave_subsection();
4574  * }
4575  *
4576  * @endcode
4577  *
4578  * Now that we have the displacement update, distribute the constraints
4579  * back to the Newton update:
4580  *
4581  * @code
4582  * constraints.distribute(newton_update);
4583  *
4584  * timer.enter_subsection("Linear solver postprocessing");
4585  * std::cout << " PP " << std::flush;
4586  *
4587  * @endcode
4588  *
4589  * The next step after solving the displacement
4590  * problem is to post-process to get the
4591  * dilatation solution from the
4592  * substitution:
4593  * @f$
4594  * d \widetilde{\mathsf{\mathbf{J}}}
4595  * = \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1} \bigl[
4596  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4597  * - \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
4598  * \bigr]
4599  * @f$
4600  *
4601  * @code
4602  * {
4603  * @endcode
4604  *
4605  * @f$
4606  * \mathsf{\mathbf{A}}_{\widetilde{p}}
4607  * =
4608  * \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
4609  * @f$
4610  *
4611  * @code
4612  * tangent_matrix.block(p_dof, u_dof)
4613  * .vmult(A.block(p_dof), newton_update.block(u_dof));
4614  * @endcode
4615  *
4616  * @f$
4617  * \mathsf{\mathbf{A}}_{\widetilde{p}}
4618  * =
4619  * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
4620  * @f$
4621  *
4622  * @code
4623  * A.block(p_dof) *= -1.0;
4624  * @endcode
4625  *
4626  * @f$
4627  * \mathsf{\mathbf{A}}_{\widetilde{p}}
4628  * =
4629  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4630  * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
4631  * @f$
4632  *
4633  * @code
4634  * A.block(p_dof) += system_rhs.block(p_dof);
4635  * @endcode
4636  *
4637  * @f$
4638  * d\mathsf{\mathbf{\widetilde{J}}}
4639  * =
4640  * \mathsf{\mathbf{K}}^{-1}_{\widetilde{p}\widetilde{J}}
4641  * [
4642  * \mathsf{\mathbf{F}}_{\widetilde{p}}
4643  * -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
4644  * ]
4645  * @f$
4646  *
4647  * @code
4648  * tangent_matrix.block(p_dof, J_dof)
4649  * .vmult(newton_update.block(J_dof), A.block(p_dof));
4650  * }
4651  *
4652  * @endcode
4653  *
4654  * we ensure here that any Dirichlet constraints
4655  * are distributed on the updated solution:
4656  *
4657  * @code
4658  * constraints.distribute(newton_update);
4659  *
4660  * @endcode
4661  *
4662  * Finally we solve for the pressure
4663  * update with the substitution:
4664  * @f$
4665  * d \widetilde{\mathsf{\mathbf{p}}}
4666  * =
4667  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4668  * \bigl[
4669  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4670  * - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4671  * d \widetilde{\mathsf{\mathbf{J}}}
4672  * \bigr]
4673  * @f$
4674  *
4675  * @code
4676  * {
4677  * @endcode
4678  *
4679  * @f$
4680  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4681  * =
4682  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4683  * d \widetilde{\mathsf{\mathbf{J}}}
4684  * @f$
4685  *
4686  * @code
4687  * tangent_matrix.block(J_dof, J_dof)
4688  * .vmult(A.block(J_dof), newton_update.block(J_dof));
4689  * @endcode
4690  *
4691  * @f$
4692  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4693  * =
4694  * -\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4695  * d \widetilde{\mathsf{\mathbf{J}}}
4696  * @f$
4697  *
4698  * @code
4699  * A.block(J_dof) *= -1.0;
4700  * @endcode
4701  *
4702  * @f$
4703  * \mathsf{\mathbf{A}}_{\widetilde{J}}
4704  * =
4705  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4706  * -
4707  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4708  * d \widetilde{\mathsf{\mathbf{J}}}
4709  * @f$
4710  *
4711  * @code
4712  * A.block(J_dof) += system_rhs.block(J_dof);
4713  * @endcode
4714  *
4715  * and finally....
4716  * @f$
4717  * d \widetilde{\mathsf{\mathbf{p}}}
4718  * =
4719  * \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
4720  * \bigl[
4721  * \mathsf{\mathbf{F}}_{\widetilde{J}}
4722  * - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
4723  * d \widetilde{\mathsf{\mathbf{J}}}
4724  * \bigr]
4725  * @f$
4726  *
4727  * @code
4728  * tangent_matrix.block(p_dof, J_dof)
4729  * .Tvmult(newton_update.block(p_dof), A.block(J_dof));
4730  * }
4731  *
4732  * @endcode
4733  *
4734  * We are now at the end, so we distribute all
4735  * constrained dofs back to the Newton
4736  * update:
4737  *
4738  * @code
4739  * constraints.distribute(newton_update);
4740  *
4741  * timer.leave_subsection();
4742  * }
4743  * else
4744  * {
4745  * std::cout << " ------ " << std::flush;
4746  *
4747  * timer.enter_subsection("Linear solver");
4748  * std::cout << " SLV " << std::flush;
4749  *
4750  * if (parameters.type_lin == "CG")
4751  * {
4752  * @endcode
4753  *
4754  * Manual condensation of the dilatation and pressure fields on
4755  * a local level, and subsequent post-processing, took quite a
4756  * bit of effort to achieve. To recap, we had to produce the
4757  * inverse matrix
4758  * @f$\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}@f$, which
4759  * was permanently written into the global tangent matrix. We then
4760  * permanently modified @f$\mathsf{\mathbf{K}}_{uu}@f$ to produce
4761  * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$. This involved the
4762  * extraction and manipulation of local sub-blocks of the tangent
4763  * matrix. After solving for the displacement, the individual
4764  * matrix-vector operations required to solve for dilatation and
4765  * pressure were carefully implemented. Contrast these many sequence
4766  * of steps to the much simpler and transparent implementation using
4767  * functionality provided by the LinearOperator class.
4768  *
4769 
4770  *
4771  * For ease of later use, we define some aliases for
4772  * blocks in the RHS vector
4773  *
4774  * @code
4775  * const Vector<double> &f_u = system_rhs.block(u_dof);
4776  * const Vector<double> &f_p = system_rhs.block(p_dof);
4777  * const Vector<double> &f_J = system_rhs.block(J_dof);
4778  *
4779  * @endcode
4780  *
4781  * ... and for blocks in the Newton update vector.
4782  *
4783  * @code
4784  * Vector<double> &d_u = newton_update.block(u_dof);
4785  * Vector<double> &d_p = newton_update.block(p_dof);
4786  * Vector<double> &d_J = newton_update.block(J_dof);
4787  *
4788  * @endcode
4789  *
4790  * We next define some linear operators for the tangent matrix
4791  * sub-blocks We will exploit the symmetry of the system, so not all
4792  * blocks are required.
4793  *
4794  * @code
4795  * const auto K_uu =
4796  * linear_operator(tangent_matrix.block(u_dof, u_dof));
4797  * const auto K_up =
4798  * linear_operator(tangent_matrix.block(u_dof, p_dof));
4799  * const auto K_pu =
4800  * linear_operator(tangent_matrix.block(p_dof, u_dof));
4801  * const auto K_Jp =
4802  * linear_operator(tangent_matrix.block(J_dof, p_dof));
4803  * const auto K_JJ =
4804  * linear_operator(tangent_matrix.block(J_dof, J_dof));
4805  *
4806  * @endcode
4807  *
4808  * We then construct a LinearOperator that represents the inverse of
4809  * (square block)
4810  * @f$\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}@f$. Since it is
4811  * diagonal (or, when a higher order ansatz it used, nearly
4812  * diagonal), a Jacobi preconditioner is suitable.
4813  *
4814  * @code
4816  * preconditioner_K_Jp_inv("jacobi");
4817  * preconditioner_K_Jp_inv.use_matrix(
4818  * tangent_matrix.block(J_dof, p_dof));
4819  * ReductionControl solver_control_K_Jp_inv(
4820  * static_cast<unsigned int>(tangent_matrix.block(J_dof, p_dof).m() *
4821  * parameters.max_iterations_lin),
4822  * 1.0e-30,
4823  * parameters.tol_lin);
4824  * SolverSelector<Vector<double>> solver_K_Jp_inv;
4825  * solver_K_Jp_inv.select("cg");
4826  * solver_K_Jp_inv.set_control(solver_control_K_Jp_inv);
4827  * const auto K_Jp_inv =
4828  * inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv);
4829  *
4830  * @endcode
4831  *
4832  * Now we can construct that transpose of
4833  * @f$\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}@f$ and a
4834  * linear operator that represents the condensed operations
4835  * @f$\overline{\mathsf{\mathbf{K}}}@f$ and
4836  * @f$\overline{\overline{\mathsf{\mathbf{K}}}}@f$ and the final
4837  * augmented matrix
4838  * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$.
4839  * Note that the schur_complement() operator could also be of use
4840  * here, but for clarity and the purpose of demonstrating the
4841  * similarities between the formulation and implementation of the
4842  * linear solution scheme, we will perform these operations
4843  * manually.
4844  *
4845  * @code
4846  * const auto K_pJ_inv = transpose_operator(K_Jp_inv);
4847  * const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv;
4848  * const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
4849  * const auto K_uu_con = K_uu + K_uu_bar_bar;
4850  *
4851  * @endcode
4852  *
4853  * Lastly, we define an operator for inverse of augmented stiffness
4854  * matrix, namely @f$\mathsf{\mathbf{K}}_{\textrm{con}}^{-1}@f$. Note
4855  * that the preconditioner for the augmented stiffness matrix is
4856  * different to the case when we use static condensation. In this
4857  * instance, the preconditioner is based on a non-modified
4858  * @f$\mathsf{\mathbf{K}}_{uu}@f$, while with the first approach we
4859  * actually modified the entries of this sub-block. However, since
4860  * @f$\mathsf{\mathbf{K}}_{\textrm{con}}@f$ and
4861  * @f$\mathsf{\mathbf{K}}_{uu}@f$ operate on the same space, it remains
4862  * adequate for this problem.
4863  *
4864  * @code
4866  * preconditioner_K_con_inv(parameters.preconditioner_type,
4867  * parameters.preconditioner_relaxation);
4868  * preconditioner_K_con_inv.use_matrix(
4869  * tangent_matrix.block(u_dof, u_dof));
4870  * ReductionControl solver_control_K_con_inv(
4871  * static_cast<unsigned int>(tangent_matrix.block(u_dof, u_dof).m() *
4872  * parameters.max_iterations_lin),
4873  * 1.0e-30,
4874  * parameters.tol_lin);
4875  * SolverSelector<Vector<double>> solver_K_con_inv;
4876  * solver_K_con_inv.select("cg");
4877  * solver_K_con_inv.set_control(solver_control_K_con_inv);
4878  * const auto K_uu_con_inv =
4879  * inverse_operator(K_uu_con,
4880  * solver_K_con_inv,
4881  * preconditioner_K_con_inv);
4882  *
4883  * @endcode
4884  *
4885  * Now we are in a position to solve for the displacement field.
4886  * We can nest the linear operations, and the result is immediately
4887  * written to the Newton update vector.
4888  * It is clear that the implementation closely mimics the derivation
4889  * stated in the introduction.
4890  *
4891  * @code
4892  * d_u =
4893  * K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p));
4894  *
4895  * timer.leave_subsection();
4896  *
4897  * @endcode
4898  *
4899  * The operations need to post-process for the dilatation and
4900  * pressure fields are just as easy to express.
4901  *
4902  * @code
4903  * timer.enter_subsection("Linear solver postprocessing");
4904  * std::cout << " PP " << std::flush;
4905  *
4906  * d_J = K_pJ_inv * (f_p - K_pu * d_u);
4907  * d_p = K_Jp_inv * (f_J - K_JJ * d_J);
4908  *
4909  * lin_it = solver_control_K_con_inv.last_step();
4910  * lin_res = solver_control_K_con_inv.last_value();
4911  * }
4912  * else if (parameters.type_lin == "Direct")
4913  * {
4914  * @endcode
4915  *
4916  * Solve the full block system with
4917  * a direct solver. As it is relatively
4918  * robust, it may be immune to problem
4919  * arising from the presence of the zero
4920  * @f$\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}@f$
4921  * block.
4922  *
4923  * @code
4924  * SparseDirectUMFPACK A_direct;
4925  * A_direct.initialize(tangent_matrix);
4926  * A_direct.vmult(newton_update, system_rhs);
4927  *
4928  * lin_it = 1;
4929  * lin_res = 0.0;
4930  *
4931  * std::cout << " -- " << std::flush;
4932  * }
4933  * else
4934  * Assert(false, ExcMessage("Linear solver type not implemented"));
4935  *
4936  * timer.leave_subsection();
4937  *
4938  * @endcode
4939  *
4940  * Finally, we again ensure here that any Dirichlet
4941  * constraints are distributed on the updated solution:
4942  *
4943  * @code
4944  * constraints.distribute(newton_update);
4945  * }
4946  *
4947  * return std::make_pair(lin_it, lin_res);
4948  * }
4949  *
4950  * @endcode
4951  *
4952  *
4953  * <a name="Solidoutput_results"></a>
4954  * <h4>Solid::output_results</h4>
4955  * Here we present how the results are written to file to be viewed
4956  * using ParaView or VisIt. The method is similar to that shown in previous
4957  * tutorials so will not be discussed in detail.
4958  *
4959  * @code
4960  * template <int dim>
4961  * void Solid<dim>::output_results() const
4962  * {
4963  * DataOut<dim> data_out;
4964  * std::vector<DataComponentInterpretation::DataComponentInterpretation>
4965  * data_component_interpretation(
4967  * data_component_interpretation.push_back(
4969  * data_component_interpretation.push_back(
4971  *
4972  * std::vector<std::string> solution_name(dim, "displacement");
4973  * solution_name.emplace_back("pressure");
4974  * solution_name.emplace_back("dilatation");
4975  *
4976  * DataOutBase::VtkFlags output_flags;
4977  * output_flags.write_higher_order_cells = true;
4978  * data_out.set_flags(output_flags);
4979  *
4980  * data_out.attach_dof_handler(dof_handler);
4981  * data_out.add_data_vector(solution_n,
4982  * solution_name,
4984  * data_component_interpretation);
4985  *
4986  * @endcode
4987  *
4988  * Since we are dealing with a large deformation problem, it would be nice
4989  * to display the result on a displaced grid! The MappingQEulerian class
4990  * linked with the DataOut class provides an interface through which this
4991  * can be achieved without physically moving the grid points in the
4992  * Triangulation object ourselves. We first need to copy the solution to
4993  * a temporary vector and then create the Eulerian mapping. We also
4994  * specify the polynomial degree to the DataOut object in order to produce
4995  * a more refined output data set when higher order polynomials are used.
4996  *
4997  * @code
4998  * Vector<double> soln(solution_n.size());
4999  * for (unsigned int i = 0; i < soln.size(); ++i)
5000  * soln(i) = solution_n(i);
5001  * MappingQEulerian<dim> q_mapping(degree, dof_handler, soln);
5002  * data_out.build_patches(q_mapping, degree);
5003  *
5004  * std::ofstream output("solution-" + std::to_string(dim) + "d-" +
5005  * std::to_string(time.get_timestep()) + ".vtu");
5006  * data_out.write_vtu(output);
5007  * }
5008  *
5009  * } // namespace Step44
5010  *
5011  *
5012  * @endcode
5013  *
5014  *
5015  * <a name="Mainfunction"></a>
5016  * <h3>Main function</h3>
5017  * Lastly we provide the main driver function which appears
5018  * no different to the other tutorials.
5019  *
5020  * @code
5021  * int main()
5022  * {
5023  * using namespace Step44;
5024  *
5025  * try
5026  * {
5027  * const unsigned int dim = 3;
5028  * Solid<dim> solid("parameters.prm");
5029  * solid.run();
5030  * }
5031  * catch (std::exception &exc)
5032  * {
5033  * std::cerr << std::endl
5034  * << std::endl
5035  * << "----------------------------------------------------"
5036  * << std::endl;
5037  * std::cerr << "Exception on processing: " << std::endl
5038  * << exc.what() << std::endl
5039  * << "Aborting!" << std::endl
5040  * << "----------------------------------------------------"
5041  * << std::endl;
5042  *
5043  * return 1;
5044  * }
5045  * catch (...)
5046  * {
5047  * std::cerr << std::endl
5048  * << std::endl
5049  * << "----------------------------------------------------"
5050  * << std::endl;
5051  * std::cerr << "Unknown exception!" << std::endl
5052  * << "Aborting!" << std::endl
5053  * << "----------------------------------------------------"
5054  * << std::endl;
5055  * return 1;
5056  * }
5057  *
5058  * return 0;
5059  * }
5060  * @endcode
5061 <a name="Results"></a><h1>Results</h1>
5062 
5063 
5064 Firstly, we present a comparison of a series of 3-d results with those
5065 in the literature (see Reese et al (2000)) to demonstrate that the program works as expected.
5066 
5067 We begin with a comparison of the convergence with mesh refinement for the @f$Q_1-DGPM_0-DGPM_0@f$ and
5068 @f$Q_2-DGPM_1-DGPM_1@f$ formulations, as summarised in the figure below.
5069 The vertical displacement of the midpoint of the upper surface of the block is used to assess convergence.
5070 Both schemes demonstrate good convergence properties for varying values of the load parameter @f$p/p_0@f$.
5071 The results agree with those in the literature.
5072 The lower-order formulation typically overestimates the displacement for low levels of refinement,
5073 while the higher-order interpolation scheme underestimates it, but be a lesser degree.
5074 This benchmark, and a series of others not shown here, give us confidence that the code is working
5075 as it should.
5076 
5077 <table align="center" class="tutorial" cellspacing="3" cellpadding="3">
5078  <tr>
5079  <td align="center">
5080  <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_convergence.png" alt="">
5081  <p align="center">
5082  Convergence of the @f$Q_1-DGPM_0-DGPM_0@f$ formulation in 3-d.
5083  </p>
5084  </td>
5085  <td align="center">
5086  <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_convergence.png" alt="">
5087  <p align="center">
5088  Convergence of the @f$Q_2-DGPM_1-DGPM_1@f$ formulation in 3-d.
5089  </p>
5090  </td>
5091  </tr>
5092 </table>
5093 
5094 
5095 A typical screen output generated by running the problem is shown below.
5096 The particular case demonstrated is that of the @f$Q_2-DGPM_1-DGPM_1@f$ formulation.
5097 It is clear that, using the Newton-Raphson method, quadratic convergence of the solution is obtained.
5098 Solution convergence is achieved within 5 Newton increments for all time-steps.
5099 The converged displacement's @f$L_2@f$-norm is several orders of magnitude less than the geometry scale.
5100 
5101 @code
5102 Grid:
5103  Reference volume: 1e-09
5105  Number of active cells: 64
5106  Number of degrees of freedom: 2699
5107  Setting up quadrature point data...
5108 
5109 Timestep 1 @ 0.1s
5110 ___________________________________________________________________________________________________________________________________________________________
5111  SOLVER STEP | LIN_IT LIN_RES RES_NORM RES_U RES_P RES_J NU_NORM NU_U NU_P NU_J
5112 ___________________________________________________________________________________________________________________________________________________________
5113  0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 786 2.118e-06 1.000e+00 1.000e+00 0.000e+00 0.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
5114  1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 552 1.031e-03 8.563e-02 8.563e-02 9.200e-13 3.929e-08 1.060e-01 3.816e-02 1.060e-01 1.060e-01
5115  2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 667 5.602e-06 2.482e-03 2.482e-03 3.373e-15 2.982e-10 2.936e-03 2.053e-04 2.936e-03 2.936e-03
5116  3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 856 6.469e-10 2.129e-06 2.129e-06 2.245e-19 1.244e-13 1.887e-06 7.289e-07 1.887e-06 1.887e-06
5117  4 ASM_R CONVERGED!
5118 ___________________________________________________________________________________________________________________________________________________________
5119 Relative errors:
5120 Displacement: 7.289e-07
5121 Force: 2.451e-10
5122 Dilatation: 1.353e-07
5123 v / V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
5124 
5125 
5126 [...]
5127 
5128 Timestep 10 @ 1.000e+00s
5129 ___________________________________________________________________________________________________________________________________________________________
5130  SOLVER STEP | LIN_IT LIN_RES RES_NORM RES_U RES_P RES_J NU_NORM NU_U NU_P NU_J
5131 ___________________________________________________________________________________________________________________________________________________________
5132  0 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 874 2.358e-06 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00 1.000e+00
5133  1 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 658 2.942e-04 1.544e-01 1.544e-01 1.208e+13 1.855e+06 6.014e-02 7.398e-02 6.014e-02 6.014e-02
5134  2 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 790 2.206e-06 2.908e-03 2.908e-03 7.302e+10 2.067e+03 2.716e-03 1.433e-03 2.716e-03 2.717e-03
5135  3 ASM_R ASM_K CST ASM_SC SLV PP UQPH | 893 2.374e-09 1.919e-06 1.919e-06 4.527e+07 4.100e+00 1.672e-06 6.842e-07 1.672e-06 1.672e-06
5136  4 ASM_R CONVERGED!
5137 ___________________________________________________________________________________________________________________________________________________________
5138 Relative errors:
5139 Displacement: 6.842e-07
5140 Force: 8.995e-10
5141 Dilatation: 1.528e-06
5142 v / V_0: 1.000e-09 / 1.000e-09 = 1.000e+00
5143 @endcode
5144 
5145 
5146 
5147 Using the Timer class, we can discern which parts of the code require the highest computational expense.
5148 For a case with a large number of degrees-of-freedom (i.e. a high level of refinement), a typical output of the Timer is given below.
5149 Much of the code in the tutorial has been developed based on the optimizations described,
5150 discussed and demonstrated in @ref step_18 "step-18" and others.
5151 With over 93% of the time being spent in the linear solver, it is obvious that it may be necessary
5152 to invest in a better solver for large three-dimensional problems.
5153 The SSOR preconditioner is not multithreaded but is effective for this class of solid problems.
5154 It may be beneficial to investigate the use of another solver such as those available through the Trilinos library.
5155 
5156 
5157 @code
5158 +---------------------------------------------+------------+------------+
5159 | Total wallclock time elapsed since start | 9.874e+02s | |
5160 | | | |
5161 | Section | no. calls | wall time | % of total |
5162 +---------------------------------+-----------+------------+------------+
5163 | Assemble system right-hand side | 53 | 1.727e+00s | 1.75e-01% |
5164 | Assemble tangent matrix | 43 | 2.707e+01s | 2.74e+00% |
5165 | Linear solver | 43 | 9.248e+02s | 9.37e+01% |
5166 | Linear solver postprocessing | 43 | 2.743e-02s | 2.78e-03% |
5167 | Perform static condensation | 43 | 1.437e+01s | 1.46e+00% |
5168 | Setup system | 1 | 3.897e-01s | 3.95e-02% |
5169 | Update QPH data | 43 | 5.770e-01s | 5.84e-02% |
5170 +---------------------------------+-----------+------------+------------+
5171 @endcode
5172 
5173 
5174 We then used ParaView to visualize the results for two cases.
5175 The first was for the coarsest grid and the lowest-order interpolation method: @f$Q_1-DGPM_0-DGPM_0@f$.
5176 The second was on a refined grid using a @f$Q_2-DGPM_1-DGPM_1@f$ formulation.
5177 The vertical component of the displacement, the pressure @f$\widetilde{p}@f$ and the dilatation @f$\widetilde{J}@f$ fields
5178 are shown below.
5179 
5180 
5181 For the first case it is clear that the coarse spatial discretization coupled with large displacements leads to a low quality solution
5182 (the loading ratio is @f$p/p_0=80@f$).
5183 Additionally, the pressure difference between elements is very large.
5184 The constant pressure field on the element means that the large pressure gradient is not captured.
5185 However, it should be noted that locking, which would be present in a standard @f$Q_1@f$ displacement formulation does not arise
5186 even in this poorly discretised case.
5187 The final vertical displacement of the tracked node on the top surface of the block is still within 12.5% of the converged solution.
5188 The pressure solution is very coarse and has large jumps between adjacent cells.
5189 It is clear that the volume nearest to the applied traction undergoes compression while the outer extents
5190 of the domain are in a state of expansion.
5191 The dilatation solution field and pressure field are clearly linked,
5192 with positive dilatation indicating regions of positive pressure and negative showing regions placed in compression.
5193 As discussed in the Introduction, a compressive pressure has a negative sign
5194 while an expansive pressure takes a positive sign.
5195 This stems from the definition of the volumetric strain energy function
5196 and is opposite to the physically realistic interpretation of pressure.
5197 
5198 
5199 <table align="center" class="tutorial" cellspacing="3" cellpadding="3">
5200  <tr>
5201  <td align="center">
5202  <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-displacement.png" alt="">
5203  <p align="center">
5204  Z-displacement solution for the 3-d problem.
5205  </p>
5206  </td>
5207  <td align="center">
5208  <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-pressure.png" alt="">
5209  <p align="center">
5210  Discontinuous piece-wise constant pressure field.
5211  </p>
5212  </td>
5213  <td align="center">
5214  <img src="https://www.dealii.org/images/steps/developer/step-44.Q1-P0_gr_1_p_ratio_80-dilatation.png" alt="">
5215  <p align="center">
5216  Discontinuous piece-wise constant dilatation field.
5217  </p>
5218  </td>
5219  </tr>
5220 </table>
5221 
5222 Combining spatial refinement and a higher-order interpolation scheme results in a high-quality solution.
5223 Three grid refinements coupled with a @f$Q_2-DGPM_1-DGPM_1@f$ formulation produces
5224 a result that clearly captures the mechanics of the problem.
5225 The deformation of the traction surface is well resolved.
5226 We can now observe the actual extent of the applied traction, with the maximum force being applied
5227 at the central point of the surface causing the largest compression.
5228 Even though very high strains are experienced in the domain,
5229 especially at the boundary of the region of applied traction,
5230 the solution remains accurate.
5231 The pressure field is captured in far greater detail than before.
5232 There is a clear distinction and transition between regions of compression and expansion,
5233 and the linear approximation of the pressure field allows a refined visualization
5234 of the pressure at the sub-element scale.
5235 It should however be noted that the pressure field remains discontinuous
5236 and could be smoothed on a continuous grid for the post-processing purposes.
5237 
5238 
5239 
5240 <table align="center" class="tutorial" cellspacing="3" cellpadding="3">
5241  <tr>
5242  <td align="center">
5243  <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-displacement.png" alt="">
5244  <p align="center">
5245  Z-displacement solution for the 3-d problem.
5246  </p>
5247  </td>
5248  <td align="center">
5249  <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-pressure.png" alt="">
5250  <p align="center">
5251  Discontinuous linear pressure field.
5252  </p>
5253  </td>
5254  <td align="center">
5255  <img src="https://www.dealii.org/images/steps/developer/step-44.Q2-P1_gr_3_p_ratio_80-dilatation.png" alt="">
5256  <p align="center">
5257  Discontinuous linear dilatation field.
5258  </p>
5259  </td>
5260  </tr>
5261 </table>
5262 
5263 This brief analysis of the results demonstrates that the three-field formulation is effective
5264 in circumventing volumetric locking for highly-incompressible media.
5265 The mixed formulation is able to accurately simulate the displacement of a
5266 near-incompressible block under compression.
5267 The command-line output indicates that the volumetric change under extreme compression resulted in
5268 less than 0.01% volume change for a Poisson's ratio of 0.4999.
5269 
5270 In terms of run-time, the @f$Q_2-DGPM_1-DGPM_1@f$ formulation tends to be more computationally expensive
5271 than the @f$Q_1-DGPM_0-DGPM_0@f$ for a similar number of degrees-of-freedom
5272 (produced by adding an extra grid refinement level for the lower-order interpolation).
5273 This is shown in the graph below for a batch of tests run consecutively on a single 4-core (8-thread) machine.
5274 The increase in computational time for the higher-order method is likely due to
5275 the increased band-width required for the higher-order elements.
5276 As previously mentioned, the use of a better solver and preconditioner may mitigate the
5277 expense of using a higher-order formulation.
5278 It was observed that for the given problem using the multithreaded Jacobi preconditioner can reduce the
5279 computational runtime by up to 72% (for the worst case being a higher-order formulation with a large number
5280 of degrees-of-freedom) in comparison to the single-thread SSOR preconditioner.
5281 However, it is the author's experience that the Jacobi method of preconditioning may not be suitable for
5282 some finite-strain problems involving alternative constitutive models.
5283 
5284 
5285 <table align="center" class="tutorial" cellspacing="3" cellpadding="3">
5286  <tr>
5287  <td align="center">
5288  <img src="https://www.dealii.org/images/steps/developer/step-44.Normalised_runtime.png" alt="">
5289  <p align="center">
5290  Runtime on a 4-core machine, normalised against the lowest grid resolution @f$Q_1-DGPM_0-DGPM_0@f$ solution that utilised a SSOR preconditioner.
5291  </p>
5292  </td>
5293  </tr>
5294 </table>
5295 
5296 
5297 Lastly, results for the displacement solution for the 2-d problem are showcased below for
5298 two different levels of grid refinement.
5299 It is clear that due to the extra constraints imposed by simulating in 2-d that the resulting
5300 displacement field, although qualitatively similar, is different to that of the 3-d case.
5301 
5302 
5303 <table align="center" class="tutorial" cellspacing="3" cellpadding="3">
5304  <tr>
5305  <td align="center">
5306  <img src="https://www.dealii.org/images/steps/developer/step-44.2d-gr_2.png" alt="">
5307  <p align="center">
5308  Y-displacement solution in 2-d for 2 global grid refinement levels.
5309  </p>
5310  </td>
5311  <td align="center">
5312  <img src="https://www.dealii.org/images/steps/developer/step-44.2d-gr_5.png" alt="">
5313  <p align="center">
5314  Y-displacement solution in 2-d for 5 global grid refinement levels.
5315  </p>
5316  </td>
5317  </tr>
5318 </table>
5319 
5320 <a name="extensions"></a>
5321 <a name="Possibilitiesforextensions"></a><h3>Possibilities for extensions</h3>
5322 
5323 
5324 There are a number of obvious extensions for this work:
5325 
5326 - Firstly, an additional constraint could be added to the free-energy
5327  function in order to enforce a high degree of incompressibility in
5328  materials. An additional Lagrange multiplier would be introduced,
5329  but this could most easily be dealt with using the principle of
5330  augmented Lagrange multipliers. This is demonstrated in <em>Simo and
5331  Taylor (1991) </em>.
5332 - The constitutive relationship used in this
5333  model is relatively basic. It may be beneficial to split the material
5334  class into two separate classes, one dealing with the volumetric
5335  response and the other the isochoric response, and produce a generic
5336  materials class (i.e. having abstract virtual functions that derived
5337  classes have to implement) that would allow for the addition of more complex
5338  material models. Such models could include other hyperelastic
5339  materials, plasticity and viscoelastic materials and others.
5340 - The program has been developed for solving problems on single-node
5341  multicore machines. With a little effort, the program could be
5342  extended to a large-scale computing environment through the use of
5343  Petsc or Trilinos, using a similar technique to that demonstrated in
5344  @ref step_40 "step-40". This would mostly involve changes to the setup, assembly,
5345  <code>PointHistory</code> and linear solver routines.
5346 - As this program assumes quasi-static equilibrium, extensions to
5347  include dynamic effects would be necessary to study problems where
5348  inertial effects are important, e.g. problems involving impact.
5349 - Load and solution limiting procedures may be necessary for highly
5350  nonlinear problems. It is possible to add a linesearch algorithm to
5351  limit the step size within a Newton increment to ensure optimum
5352  convergence. It may also be necessary to use a load limiting method,
5353  such as the Riks method, to solve unstable problems involving
5354  geometric instability such as buckling and snap-through.
5355 - Many physical problems involve contact. It is possible to include
5356  the effect of frictional or frictionless contact between objects
5357  into this program. This would involve the addition of an extra term
5358  in the free-energy functional and therefore an addition to the
5359  assembly routine. One would also need to manage the contact problem
5360  (detection and stress calculations) itself. An alternative to
5361  additional penalty terms in the free-energy functional would be to
5362  use active set methods such as the one used in @ref step_41 "step-41".
5363 - The complete condensation procedure using LinearOperators has been
5364  coded into the linear solver routine. This could also have been
5365  achieved through the application of the schur_complement()
5366  operator to condense out one or more of the fields in a more
5367  automated manner.
5368 - Finally, adaptive mesh refinement, as demonstrated in @ref step_6 "step-6" and
5369  @ref step_18 "step-18", could provide additional solution accuracy.
5370  *
5371  *
5372 <a name="PlainProg"></a>
5373 <h1> The plain program</h1>
5374 @include "step-44.cc"
5375 */
Physics::Elasticity::Kinematics::F
Tensor< 2, dim, Number > F(const Tensor< 2, dim, Number > &Grad_u)
LAPACKSupport::diagonal
@ diagonal
Matrix is diagonal.
Definition: lapack_support.h:121
TimerOutput::enter_subsection
void enter_subsection(const std::string &section_name)
Definition: timer.cc:403
FEValuesExtractors
Definition: fe_values_extractors.h:82
DataOutInterface::write_vtu
void write_vtu(std::ostream &out) const
Definition: data_out_base.cc:6864
SolverSelector
Definition: solver_selector.h:91
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria, const Mapping< dim, spacedim > &mapping=(StaticMappingQ1< dim, spacedim >::mapping))
Definition: grid_tools.cc:133
SparseDirectUMFPACK::initialize
void initialize(const SparsityPattern &sparsity_pattern)
Definition: sparse_direct.cc:49
SymmetricTensor< 2, dim >
BlockVector< double >
DataComponentInterpretation::component_is_scalar
@ component_is_scalar
Definition: data_component_interpretation.h:55
WorkStream
Definition: work_stream.h:157
Triangulation< dim >
FEValuesExtractors::Scalar
Definition: fe_values_extractors.h:95
Patterns::Tools::to_string
std::string to_string(const T &t)
Definition: patterns.h:2360
PreconditionSelector
Definition: precondition_selector.h:103
Physics::Elasticity::Kinematics::C
SymmetricTensor< 2, dim, Number > C(const Tensor< 2, dim, Number > &F)
DoFHandler::n_components
unsigned int n_components(const DoFHandler< dim, spacedim > &dh)
ComponentMask
Definition: component_mask.h:83
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >(0)), const bool project_to_boundary_first=false)
ReductionControl
Definition: solver_control.h:427
types::boundary_id
unsigned int boundary_id
Definition: types.h:129
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
second
Point< 2 > second
Definition: grid_out.cc:4353
LinearOperator::inverse_operator
LinearOperator< Domain, Range, Payload > inverse_operator(const LinearOperator< Range, Domain, Payload > &op, Solver &solver, const Preconditioner &preconditioner)
Definition: linear_operator.h:693
DataOut::build_patches
virtual void build_patches(const unsigned int n_subdivisions=0)
Definition: data_out.cc:1071
Threads::new_task
Task< RT > new_task(const std::function< RT()> &function)
Definition: thread_management.h:1647
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
LinearAlgebra::CUDAWrappers::kernel::set
__global__ void set(Number *val, const Number s, const size_type N)
LinearOperator::schur_complement
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)
Definition: schur_complement.h:249
OpenCASCADE::point
Point< spacedim > point(const gp_Pnt &p, const double tolerance=1e-10)
Definition: utilities.cc:188
FEValues< dim >
Timer
Definition: timer.h:119
Utilities::CUDA::free
void free(T *&pointer)
Definition: cuda.h:96
WorkStream::run
void run(const std::vector< std::vector< Iterator >> &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition: work_stream.h:1185
level
unsigned int level
Definition: grid_out.cc:4355
TimerOutput
Definition: timer.h:546
LinearOperator::linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const Matrix &matrix)
Definition: linear_operator.h:1373
LAPACKSupport::one
static const types::blas_int one
Definition: lapack_support.h:183
DataOutBase::none
@ none
Definition: data_out_base.h:1557
LAPACKSupport::symmetric
@ symmetric
Matrix is symmetric.
Definition: lapack_support.h:115
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandlerType< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask=ComponentMask())
SolverSelector::select
void select(const std::string &name)
Definition: solver_selector.h:261
Algorithms::Events::initial
const Event initial
Definition: event.cc:65
StandardExceptions::ExcMessage
static ::ExceptionBase & ExcMessage(std::string arg1)
TrilinosWrappers::internal::begin
VectorType::value_type * begin(VectorType &V)
Definition: trilinos_sparse_matrix.cc:51
TimerOutput::leave_subsection
void leave_subsection(const std::string &section_name="")
Definition: timer.cc:445
Tensor< 2, dim >
ComponentSelectFunction
Definition: function.h:560
LocalIntegrators::Divergence::norm
double norm(const FEValuesBase< dim > &fe, const ArrayView< const std::vector< Tensor< 1, dim >>> &Du)
Definition: divergence.h:548
GridTools::scale
void scale(const double scaling_factor, Triangulation< dim, spacedim > &triangulation)
Definition: grid_tools.cc:837
LAPACKSupport::matrix
@ matrix
Contents is actually a matrix.
Definition: lapack_support.h:60
LAPACKSupport::general
@ general
No special properties.
Definition: lapack_support.h:113
TrilinosWrappers::internal::end
VectorType::value_type * end(VectorType &V)
Definition: trilinos_sparse_matrix.cc:65
SolverSelector::set_control
void set_control(SolverControl &ctrl)
Definition: solver_selector.h:314
FiniteElement< dim >
UpdateFlags
UpdateFlags
Definition: fe_update_flags.h:66
DataOutBase::VtkFlags
Definition: data_out_base.h:1095
QGauss
Definition: quadrature_lib.h:40
LAPACKSupport::A
static const char A
Definition: lapack_support.h:155
internal::dummy
const types::global_dof_index * dummy()
Definition: dof_handler.cc:59
DataOut_DoFData::attach_dof_handler
void attach_dof_handler(const DoFHandlerType &)
StandardExceptions::ExcInternalError
static ::ExceptionBase & ExcInternalError()
AffineConstraints< double >
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
LinearOperator::transpose_operator
LinearOperator< Domain, Range, Payload > transpose_operator(const LinearOperator< Range, Domain, Payload > &op)
Definition: linear_operator.h:651
internal::assemble
void assemble(const MeshWorker::DoFInfoBox< dim, DOFINFO > &dinfo, A *assembler)
Definition: loop.h:71
DataOutInterface::set_flags
void set_flags(const FlagType &flags)
Definition: data_out_base.cc:7837
LAPACKSupport::zero
static const types::blas_int zero
Definition: lapack_support.h:179
DataOutBase::VtkFlags::write_higher_order_cells
bool write_higher_order_cells
Definition: data_out_base.h:1178
Functions::ZeroFunction
Definition: function.h:513
MappingQEulerian
Definition: mapping_q_eulerian.h:95
Differentiation::SD::sign
Expression sign(const Expression &x)
Definition: symengine_math.cc:280
SparseDirectUMFPACK
Definition: sparse_direct.h:90
FullMatrix< double >
triangulation
const typename ::parallel::distributed::Triangulation< dim, spacedim > * triangulation
Definition: p4est_wrappers.cc:69
FEFaceValues< dim >
MeshWorker::loop
void loop(ITERATOR begin, typename identity< ITERATOR >::type end, DOFINFO &dinfo, INFOBOX &info, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &cell_worker, const std::function< void(DOFINFO &, typename INFOBOX::CellInfo &)> &boundary_worker, const std::function< void(DOFINFO &, DOFINFO &, typename INFOBOX::CellInfo &, typename INFOBOX::CellInfo &)> &face_worker, ASSEMBLER &assembler, const LoopControl &lctrl=LoopControl())
Definition: loop.h:443
first
Point< 2 > first
Definition: grid_out.cc:4352
DataOut< dim >
Vector< double >
AffineConstraints::close
void close()
center
Point< 3 > center
Definition: data_out_base.cc:169
LinearOperator
Definition: linear_operator.h:169
DerivativeForm::transpose
DerivativeForm< 1, spacedim, dim, Number > transpose(const DerivativeForm< 1, dim, spacedim, Number > &DF)
Definition: derivative_form.h:470
internal::VectorOperations::copy
void copy(const T *begin, const T *end, U *dest)
Definition: vector_operations_internal.h:67
DataComponentInterpretation::component_is_part_of_vector
@ component_is_part_of_vector
Definition: data_component_interpretation.h:61
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