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