1 /* ---------------------------------------------------------------------
2 *
3 * Copyright (C) 2010 - 2020 by the deal.II authors and
4 * & Jean-Paul Pelteret and Andrew McBride
5 *
6 * This file is part of the deal.II library.
7 *
8 * The deal.II library is free software; you can use it, redistribute
9 * it, and/or modify it under the terms of the GNU Lesser General
10 * Public License as published by the Free Software Foundation; either
11 * version 2.1 of the License, or (at your option) any later version.
12 * The full text of the license can be found in the file LICENSE.md at
13 * the top level directory of deal.II.
14 *
15 * ---------------------------------------------------------------------
16
17 *
18 * Authors: Jean-Paul Pelteret, University of Cape Town,
19 * Andrew McBride, University of Erlangen-Nuremberg, 2010
20 */
21
22
23 // We start by including all the necessary deal.II header files and some C++
24 // related ones. They have been discussed in detail in previous tutorial
25 // programs, so you need only refer to past tutorials for details.
26 #include <deal.II/base/function.h>
27 #include <deal.II/base/parameter_handler.h>
28 #include <deal.II/base/point.h>
29 #include <deal.II/base/quadrature_lib.h>
30 #include <deal.II/base/symmetric_tensor.h>
31 #include <deal.II/base/tensor.h>
32 #include <deal.II/base/timer.h>
33 #include <deal.II/base/work_stream.h>
34 #include <deal.II/dofs/dof_renumbering.h>
35 #include <deal.II/dofs/dof_tools.h>
36
37 // This header gives us the functionality to store
38 // data at quadrature points
39 #include <deal.II/base/quadrature_point_data.h>
40
41 #include <deal.II/grid/grid_generator.h>
42 #include <deal.II/grid/grid_tools.h>
43 #include <deal.II/grid/grid_in.h>
44 #include <deal.II/grid/tria.h>
45
46 #include <deal.II/fe/fe_dgp_monomial.h>
47 #include <deal.II/fe/fe_q.h>
48 #include <deal.II/fe/fe_system.h>
49 #include <deal.II/fe/fe_tools.h>
50 #include <deal.II/fe/fe_values.h>
51 #include <deal.II/fe/mapping_q_eulerian.h>
52
53 #include <deal.II/lac/block_sparse_matrix.h>
54 #include <deal.II/lac/block_vector.h>
55 #include <deal.II/lac/dynamic_sparsity_pattern.h>
56 #include <deal.II/lac/full_matrix.h>
57 #include <deal.II/lac/precondition_selector.h>
58 #include <deal.II/lac/solver_cg.h>
59 #include <deal.II/lac/solver_selector.h>
60 #include <deal.II/lac/sparse_direct.h>
61 #include <deal.II/lac/affine_constraints.h>
62
63 // Here are the headers necessary to use the LinearOperator class.
64 // These are also all conveniently packaged into a single
65 // header file, namely <deal.II/lac/linear_operator_tools.h>
66 // but we list those specifically required here for the sake
67 // of transparency.
68 #include <deal.II/lac/linear_operator.h>
69 #include <deal.II/lac/packaged_operation.h>
70
71 #include <deal.II/numerics/data_out.h>
72 #include <deal.II/numerics/vector_tools.h>
73
74 // Defined in these two headers are some operations that are pertinent to
75 // finite strain elasticity. The first will help us compute some kinematic
76 // quantities, and the second provides some stanard tensor definitions.
77 #include <deal.II/physics/elasticity/kinematics.h>
78 #include <deal.II/physics/elasticity/standard_tensors.h>
79
80 #include <iostream>
81 #include <fstream>
82
83
84 // We then stick everything that relates to this tutorial program into a
85 // namespace of its own, and import all the deal.II function and class names
86 // into it:
87 namespace Step44
88 {
89 using namespace dealii;
90
91 // @sect3{Run-time parameters}
92 //
93 // There are several parameters that can be set in the code so we set up a
94 // ParameterHandler object to read in the choices at run-time.
95 namespace Parameters
96 {
97 // @sect4{Finite Element system}
98
99 // As mentioned in the introduction, a different order interpolation should
100 // be used for the displacement $\mathbf{u}$ than for the pressure
101 // $\widetilde{p}$ and the dilatation $\widetilde{J}$. Choosing
102 // $\widetilde{p}$ and $\widetilde{J}$ as discontinuous (constant) functions
103 // at the element level leads to the mean-dilatation method. The
104 // discontinuous approximation allows $\widetilde{p}$ and $\widetilde{J}$ to
105 // be condensed out and a classical displacement based method is recovered.
106 // Here we specify the polynomial order used to approximate the solution.
107 // The quadrature order should be adjusted accordingly.
108 struct FESystem
109 {
110 unsigned int poly_degree;
111 unsigned int quad_order;
112
113 static void declare_parameters(ParameterHandler &prm);
114
115 void parse_parameters(ParameterHandler &prm);
116 };
117
118
declare_parameters(ParameterHandler & prm)119 void FESystem::declare_parameters(ParameterHandler &prm)
120 {
121 prm.enter_subsection("Finite element system");
122 {
123 prm.declare_entry("Polynomial degree",
124 "2",
125 Patterns::Integer(0),
126 "Displacement system polynomial order");
127
128 prm.declare_entry("Quadrature order",
129 "3",
130 Patterns::Integer(0),
131 "Gauss quadrature order");
132 }
133 prm.leave_subsection();
134 }
135
parse_parameters(ParameterHandler & prm)136 void FESystem::parse_parameters(ParameterHandler &prm)
137 {
138 prm.enter_subsection("Finite element system");
139 {
140 poly_degree = prm.get_integer("Polynomial degree");
141 quad_order = prm.get_integer("Quadrature order");
142 }
143 prm.leave_subsection();
144 }
145
146 // @sect4{Geometry}
147
148 // Make adjustments to the problem geometry and the applied load. Since the
149 // problem modelled here is quite specific, the load scale can be altered to
150 // specific values to compare with the results given in the literature.
151 struct Geometry
152 {
153 unsigned int global_refinement;
154 double scale;
155 double p_p0;
156
157 static void declare_parameters(ParameterHandler &prm);
158
159 void parse_parameters(ParameterHandler &prm);
160 };
161
declare_parameters(ParameterHandler & prm)162 void Geometry::declare_parameters(ParameterHandler &prm)
163 {
164 prm.enter_subsection("Geometry");
165 {
166 prm.declare_entry("Global refinement",
167 "2",
168 Patterns::Integer(0),
169 "Global refinement level");
170
171 prm.declare_entry("Grid scale",
172 "1e-3",
173 Patterns::Double(0.0),
174 "Global grid scaling factor");
175
176 prm.declare_entry("Pressure ratio p/p0",
177 "100",
178 Patterns::Selection("20|40|60|80|100"),
179 "Ratio of applied pressure to reference pressure");
180 }
181 prm.leave_subsection();
182 }
183
parse_parameters(ParameterHandler & prm)184 void Geometry::parse_parameters(ParameterHandler &prm)
185 {
186 prm.enter_subsection("Geometry");
187 {
188 global_refinement = prm.get_integer("Global refinement");
189 scale = prm.get_double("Grid scale");
190 p_p0 = prm.get_double("Pressure ratio p/p0");
191 }
192 prm.leave_subsection();
193 }
194
195 // @sect4{Materials}
196
197 // We also need the shear modulus $ \mu $ and Poisson ration $ \nu $ for the
198 // neo-Hookean material.
199 struct Materials
200 {
201 double nu;
202 double mu;
203
204 static void declare_parameters(ParameterHandler &prm);
205
206 void parse_parameters(ParameterHandler &prm);
207 };
208
declare_parameters(ParameterHandler & prm)209 void Materials::declare_parameters(ParameterHandler &prm)
210 {
211 prm.enter_subsection("Material properties");
212 {
213 prm.declare_entry("Poisson's ratio",
214 "0.4999",
215 Patterns::Double(-1.0, 0.5),
216 "Poisson's ratio");
217
218 prm.declare_entry("Shear modulus",
219 "80.194e6",
220 Patterns::Double(),
221 "Shear modulus");
222 }
223 prm.leave_subsection();
224 }
225
parse_parameters(ParameterHandler & prm)226 void Materials::parse_parameters(ParameterHandler &prm)
227 {
228 prm.enter_subsection("Material properties");
229 {
230 nu = prm.get_double("Poisson's ratio");
231 mu = prm.get_double("Shear modulus");
232 }
233 prm.leave_subsection();
234 }
235
236 // @sect4{Linear solver}
237
238 // Next, we choose both solver and preconditioner settings. The use of an
239 // effective preconditioner is critical to ensure convergence when a large
240 // nonlinear motion occurs within a Newton increment.
241 struct LinearSolver
242 {
243 std::string type_lin;
244 double tol_lin;
245 double max_iterations_lin;
246 bool use_static_condensation;
247 std::string preconditioner_type;
248 double preconditioner_relaxation;
249
250 static void declare_parameters(ParameterHandler &prm);
251
252 void parse_parameters(ParameterHandler &prm);
253 };
254
declare_parameters(ParameterHandler & prm)255 void LinearSolver::declare_parameters(ParameterHandler &prm)
256 {
257 prm.enter_subsection("Linear solver");
258 {
259 prm.declare_entry("Solver type",
260 "CG",
261 Patterns::Selection("CG|Direct"),
262 "Type of solver used to solve the linear system");
263
264 prm.declare_entry("Residual",
265 "1e-6",
266 Patterns::Double(0.0),
267 "Linear solver residual (scaled by residual norm)");
268
269 prm.declare_entry(
270 "Max iteration multiplier",
271 "1",
272 Patterns::Double(0.0),
273 "Linear solver iterations (multiples of the system matrix size)");
274
275 prm.declare_entry("Use static condensation",
276 "true",
277 Patterns::Bool(),
278 "Solve the full block system or a reduced problem");
279
280 prm.declare_entry("Preconditioner type",
281 "ssor",
282 Patterns::Selection("jacobi|ssor"),
283 "Type of preconditioner");
284
285 prm.declare_entry("Preconditioner relaxation",
286 "0.65",
287 Patterns::Double(0.0),
288 "Preconditioner relaxation value");
289 }
290 prm.leave_subsection();
291 }
292
parse_parameters(ParameterHandler & prm)293 void LinearSolver::parse_parameters(ParameterHandler &prm)
294 {
295 prm.enter_subsection("Linear solver");
296 {
297 type_lin = prm.get("Solver type");
298 tol_lin = prm.get_double("Residual");
299 max_iterations_lin = prm.get_double("Max iteration multiplier");
300 use_static_condensation = prm.get_bool("Use static condensation");
301 preconditioner_type = prm.get("Preconditioner type");
302 preconditioner_relaxation = prm.get_double("Preconditioner relaxation");
303 }
304 prm.leave_subsection();
305 }
306
307 // @sect4{Nonlinear solver}
308
309 // A Newton-Raphson scheme is used to solve the nonlinear system of
310 // governing equations. We now define the tolerances and the maximum number
311 // of iterations for the Newton-Raphson nonlinear solver.
312 struct NonlinearSolver
313 {
314 unsigned int max_iterations_NR;
315 double tol_f;
316 double tol_u;
317
318 static void declare_parameters(ParameterHandler &prm);
319
320 void parse_parameters(ParameterHandler &prm);
321 };
322
declare_parameters(ParameterHandler & prm)323 void NonlinearSolver::declare_parameters(ParameterHandler &prm)
324 {
325 prm.enter_subsection("Nonlinear solver");
326 {
327 prm.declare_entry("Max iterations Newton-Raphson",
328 "10",
329 Patterns::Integer(0),
330 "Number of Newton-Raphson iterations allowed");
331
332 prm.declare_entry("Tolerance force",
333 "1.0e-9",
334 Patterns::Double(0.0),
335 "Force residual tolerance");
336
337 prm.declare_entry("Tolerance displacement",
338 "1.0e-6",
339 Patterns::Double(0.0),
340 "Displacement error tolerance");
341 }
342 prm.leave_subsection();
343 }
344
parse_parameters(ParameterHandler & prm)345 void NonlinearSolver::parse_parameters(ParameterHandler &prm)
346 {
347 prm.enter_subsection("Nonlinear solver");
348 {
349 max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
350 tol_f = prm.get_double("Tolerance force");
351 tol_u = prm.get_double("Tolerance displacement");
352 }
353 prm.leave_subsection();
354 }
355
356 // @sect4{Time}
357
358 // Set the timestep size $ \varDelta t $ and the simulation end-time.
359 struct Time
360 {
361 double delta_t;
362 double end_time;
363
364 static void declare_parameters(ParameterHandler &prm);
365
366 void parse_parameters(ParameterHandler &prm);
367 };
368
declare_parameters(ParameterHandler & prm)369 void Time::declare_parameters(ParameterHandler &prm)
370 {
371 prm.enter_subsection("Time");
372 {
373 prm.declare_entry("End time", "1", Patterns::Double(), "End time");
374
375 prm.declare_entry("Time step size",
376 "0.1",
377 Patterns::Double(),
378 "Time step size");
379 }
380 prm.leave_subsection();
381 }
382
parse_parameters(ParameterHandler & prm)383 void Time::parse_parameters(ParameterHandler &prm)
384 {
385 prm.enter_subsection("Time");
386 {
387 end_time = prm.get_double("End time");
388 delta_t = prm.get_double("Time step size");
389 }
390 prm.leave_subsection();
391 }
392
393 // @sect4{All parameters}
394
395 // Finally we consolidate all of the above structures into a single
396 // container that holds all of our run-time selections.
397 struct AllParameters : public FESystem,
398 public Geometry,
399 public Materials,
400 public LinearSolver,
401 public NonlinearSolver,
402 public Time
403
404 {
405 AllParameters(const std::string &input_file);
406
407 static void declare_parameters(ParameterHandler &prm);
408
409 void parse_parameters(ParameterHandler &prm);
410 };
411
AllParameters(const std::string & input_file)412 AllParameters::AllParameters(const std::string &input_file)
413 {
414 ParameterHandler prm;
415 declare_parameters(prm);
416 prm.parse_input(input_file);
417 parse_parameters(prm);
418 }
419
declare_parameters(ParameterHandler & prm)420 void AllParameters::declare_parameters(ParameterHandler &prm)
421 {
422 FESystem::declare_parameters(prm);
423 Geometry::declare_parameters(prm);
424 Materials::declare_parameters(prm);
425 LinearSolver::declare_parameters(prm);
426 NonlinearSolver::declare_parameters(prm);
427 Time::declare_parameters(prm);
428 }
429
parse_parameters(ParameterHandler & prm)430 void AllParameters::parse_parameters(ParameterHandler &prm)
431 {
432 FESystem::parse_parameters(prm);
433 Geometry::parse_parameters(prm);
434 Materials::parse_parameters(prm);
435 LinearSolver::parse_parameters(prm);
436 NonlinearSolver::parse_parameters(prm);
437 Time::parse_parameters(prm);
438 }
439 } // namespace Parameters
440
441 // @sect3{Time class}
442
443 // A simple class to store time data. Its functioning is transparent so no
444 // discussion is necessary. For simplicity we assume a constant time step
445 // size.
446 class Time
447 {
448 public:
Time(const double time_end,const double delta_t)449 Time(const double time_end, const double delta_t)
450 : timestep(0)
451 , time_current(0.0)
452 , time_end(time_end)
453 , delta_t(delta_t)
454 {}
455
456 virtual ~Time() = default;
457
current() const458 double current() const
459 {
460 return time_current;
461 }
end() const462 double end() const
463 {
464 return time_end;
465 }
get_delta_t() const466 double get_delta_t() const
467 {
468 return delta_t;
469 }
get_timestep() const470 unsigned int get_timestep() const
471 {
472 return timestep;
473 }
increment()474 void increment()
475 {
476 time_current += delta_t;
477 ++timestep;
478 }
479
480 private:
481 unsigned int timestep;
482 double time_current;
483 const double time_end;
484 const double delta_t;
485 };
486
487 // @sect3{Compressible neo-Hookean material within a three-field formulation}
488
489 // As discussed in the Introduction, Neo-Hookean materials are a type of
490 // hyperelastic materials. The entire domain is assumed to be composed of a
491 // compressible neo-Hookean material. This class defines the behaviour of
492 // this material within a three-field formulation. Compressible neo-Hookean
493 // materials can be described by a strain-energy function (SEF) $ \Psi =
494 // \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(\widetilde{J})
495 // $.
496 //
497 // The isochoric response is given by $
498 // \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] $
499 // where $ c_{1} = \frac{\mu}{2} $ and $\overline{I}_{1}$ is the first
500 // invariant of the left- or right-isochoric Cauchy-Green deformation tensors.
501 // That is $\overline{I}_1 \dealcoloneq \textrm{tr}(\overline{\mathbf{b}})$.
502 // In this example the SEF that governs the volumetric response is defined as
503 // $ \Psi_{\text{vol}}(\widetilde{J}) = \kappa \frac{1}{4} [ \widetilde{J}^2 -
504 // 1 - 2\textrm{ln}\; \widetilde{J} ]$, where $\kappa \dealcoloneq \lambda +
505 // 2/3 \mu$ is the <a href="http://en.wikipedia.org/wiki/Bulk_modulus">bulk
506 // modulus</a> and $\lambda$ is <a
507 // href="http://en.wikipedia.org/wiki/Lam%C3%A9_parameters">Lamé's
508 // first parameter</a>.
509 //
510 // The following class will be used to characterize the material we work with,
511 // and provides a central point that one would need to modify if one were to
512 // implement a different material model. For it to work, we will store one
513 // object of this type per quadrature point, and in each of these objects
514 // store the current state (characterized by the values or measures of the
515 // three fields) so that we can compute the elastic coefficients linearized
516 // around the current state.
517 template <int dim>
518 class Material_Compressible_Neo_Hook_Three_Field
519 {
520 public:
Material_Compressible_Neo_Hook_Three_Field(const double mu,const double nu)521 Material_Compressible_Neo_Hook_Three_Field(const double mu, const double nu)
522 : kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu)))
523 , c_1(mu / 2.0)
524 , det_F(1.0)
525 , p_tilde(0.0)
526 , J_tilde(1.0)
527 , b_bar(Physics::Elasticity::StandardTensors<dim>::I)
528 {
529 Assert(kappa > 0, ExcInternalError());
530 }
531
532 // We update the material model with various deformation dependent data
533 // based on $F$ and the pressure $\widetilde{p}$ and dilatation
534 // $\widetilde{J}$, and at the end of the function include a physical
535 // check for internal consistency:
update_material_data(const Tensor<2,dim> & F,const double p_tilde_in,const double J_tilde_in)536 void update_material_data(const Tensor<2, dim> &F,
537 const double p_tilde_in,
538 const double J_tilde_in)
539 {
540 det_F = determinant(F);
541 const Tensor<2, dim> F_bar = Physics::Elasticity::Kinematics::F_iso(F);
542 b_bar = Physics::Elasticity::Kinematics::b(F_bar);
543 p_tilde = p_tilde_in;
544 J_tilde = J_tilde_in;
545
546 Assert(det_F > 0, ExcInternalError());
547 }
548
549 // The second function determines the Kirchhoff stress $\boldsymbol{\tau}
550 // = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}$
get_tau()551 SymmetricTensor<2, dim> get_tau()
552 {
553 return get_tau_iso() + get_tau_vol();
554 }
555
556 // The fourth-order elasticity tensor in the spatial setting
557 // $\mathfrak{c}$ is calculated from the SEF $\Psi$ as $ J
558 // \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}$
559 // where $ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial
560 // \mathbf{C} \partial \mathbf{C}}$
get_Jc() const561 SymmetricTensor<4, dim> get_Jc() const
562 {
563 return get_Jc_vol() + get_Jc_iso();
564 }
565
566 // Derivative of the volumetric free energy with respect to
567 // $\widetilde{J}$ return $\frac{\partial
568 // \Psi_{\text{vol}}(\widetilde{J})}{\partial \widetilde{J}}$
get_dPsi_vol_dJ() const569 double get_dPsi_vol_dJ() const
570 {
571 return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde);
572 }
573
574 // Second derivative of the volumetric free energy wrt $\widetilde{J}$. We
575 // need the following computation explicitly in the tangent so we make it
576 // public. We calculate $\frac{\partial^2
577 // \Psi_{\textrm{vol}}(\widetilde{J})}{\partial \widetilde{J} \partial
578 // \widetilde{J}}$
get_d2Psi_vol_dJ2() const579 double get_d2Psi_vol_dJ2() const
580 {
581 return ((kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde)));
582 }
583
584 // The next few functions return various data that we choose to store with
585 // the material:
get_det_F() const586 double get_det_F() const
587 {
588 return det_F;
589 }
590
get_p_tilde() const591 double get_p_tilde() const
592 {
593 return p_tilde;
594 }
595
get_J_tilde() const596 double get_J_tilde() const
597 {
598 return J_tilde;
599 }
600
601 protected:
602 // Define constitutive model parameters $\kappa$ (bulk modulus) and the
603 // neo-Hookean model parameter $c_1$:
604 const double kappa;
605 const double c_1;
606
607 // Model specific data that is convenient to store with the material:
608 double det_F;
609 double p_tilde;
610 double J_tilde;
611 SymmetricTensor<2, dim> b_bar;
612
613 // The following functions are used internally in determining the result
614 // of some of the public functions above. The first one determines the
615 // volumetric Kirchhoff stress $\boldsymbol{\tau}_{\textrm{vol}}$:
get_tau_vol() const616 SymmetricTensor<2, dim> get_tau_vol() const
617 {
618 return p_tilde * det_F * Physics::Elasticity::StandardTensors<dim>::I;
619 }
620
621 // Next, determine the isochoric Kirchhoff stress
622 // $\boldsymbol{\tau}_{\textrm{iso}} =
623 // \mathcal{P}:\overline{\boldsymbol{\tau}}$:
get_tau_iso() const624 SymmetricTensor<2, dim> get_tau_iso() const
625 {
626 return Physics::Elasticity::StandardTensors<dim>::dev_P * get_tau_bar();
627 }
628
629 // Then, determine the fictitious Kirchhoff stress
630 // $\overline{\boldsymbol{\tau}}$:
get_tau_bar() const631 SymmetricTensor<2, dim> get_tau_bar() const
632 {
633 return 2.0 * c_1 * b_bar;
634 }
635
636 // Calculate the volumetric part of the tangent $J
637 // \mathfrak{c}_\textrm{vol}$:
get_Jc_vol() const638 SymmetricTensor<4, dim> get_Jc_vol() const
639 {
640 return p_tilde * det_F *
641 (Physics::Elasticity::StandardTensors<dim>::IxI -
642 (2.0 * Physics::Elasticity::StandardTensors<dim>::S));
643 }
644
645 // Calculate the isochoric part of the tangent $J
646 // \mathfrak{c}_\textrm{iso}$:
get_Jc_iso() const647 SymmetricTensor<4, dim> get_Jc_iso() const
648 {
649 const SymmetricTensor<2, dim> tau_bar = get_tau_bar();
650 const SymmetricTensor<2, dim> tau_iso = get_tau_iso();
651 const SymmetricTensor<4, dim> tau_iso_x_I =
652 outer_product(tau_iso, Physics::Elasticity::StandardTensors<dim>::I);
653 const SymmetricTensor<4, dim> I_x_tau_iso =
654 outer_product(Physics::Elasticity::StandardTensors<dim>::I, tau_iso);
655 const SymmetricTensor<4, dim> c_bar = get_c_bar();
656
657 return (2.0 / dim) * trace(tau_bar) *
658 Physics::Elasticity::StandardTensors<dim>::dev_P -
659 (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) +
660 Physics::Elasticity::StandardTensors<dim>::dev_P * c_bar *
661 Physics::Elasticity::StandardTensors<dim>::dev_P;
662 }
663
664 // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$.
665 // For the material model chosen this is simply zero:
get_c_bar() const666 SymmetricTensor<4, dim> get_c_bar() const
667 {
668 return SymmetricTensor<4, dim>();
669 }
670 };
671
672 // @sect3{Quadrature point history}
673
674 // As seen in step-18, the <code> PointHistory </code> class offers a method
675 // for storing data at the quadrature points. Here each quadrature point
676 // holds a pointer to a material description. Thus, different material models
677 // can be used in different regions of the domain. Among other data, we
678 // choose to store the Kirchhoff stress $\boldsymbol{\tau}$ and the tangent
679 // $J\mathfrak{c}$ for the quadrature points.
680 template <int dim>
681 class PointHistory
682 {
683 public:
PointHistory()684 PointHistory()
685 : F_inv(Physics::Elasticity::StandardTensors<dim>::I)
686 , tau(SymmetricTensor<2, dim>())
687 , d2Psi_vol_dJ2(0.0)
688 , dPsi_vol_dJ(0.0)
689 , Jc(SymmetricTensor<4, dim>())
690 {}
691
692 virtual ~PointHistory() = default;
693
694 // The first function is used to create a material object and to
695 // initialize all tensors correctly: The second one updates the stored
696 // values and stresses based on the current deformation measure
697 // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$, pressure $\widetilde{p}$ and
698 // dilation $\widetilde{J}$ field values.
setup_lqp(const Parameters::AllParameters & parameters)699 void setup_lqp(const Parameters::AllParameters ¶meters)
700 {
701 material =
702 std::make_shared<Material_Compressible_Neo_Hook_Three_Field<dim>>(
703 parameters.mu, parameters.nu);
704 update_values(Tensor<2, dim>(), 0.0, 1.0);
705 }
706
707 // To this end, we calculate the deformation gradient $\mathbf{F}$ from
708 // the displacement gradient $\textrm{Grad}\ \mathbf{u}$, i.e.
709 // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + \textrm{Grad}\ \mathbf{u}$ and
710 // then let the material model associated with this quadrature point
711 // update itself. When computing the deformation gradient, we have to take
712 // care with which data types we compare the sum $\mathbf{I} +
713 // \textrm{Grad}\ \mathbf{u}$: Since $I$ has data type SymmetricTensor,
714 // just writing <code>I + Grad_u_n</code> would convert the second
715 // argument to a symmetric tensor, perform the sum, and then cast the
716 // result to a Tensor (i.e., the type of a possibly nonsymmetric
717 // tensor). However, since <code>Grad_u_n</code> is nonsymmetric in
718 // general, the conversion to SymmetricTensor will fail. We can avoid this
719 // back and forth by converting $I$ to Tensor first, and then performing
720 // the addition as between nonsymmetric tensors:
update_values(const Tensor<2,dim> & Grad_u_n,const double p_tilde,const double J_tilde)721 void update_values(const Tensor<2, dim> &Grad_u_n,
722 const double p_tilde,
723 const double J_tilde)
724 {
725 const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(Grad_u_n);
726 material->update_material_data(F, p_tilde, J_tilde);
727
728 // The material has been updated so we now calculate the Kirchhoff
729 // stress $\mathbf{\tau}$, the tangent $J\mathfrak{c}$ and the first and
730 // second derivatives of the volumetric free energy.
731 //
732 // We also store the inverse of the deformation gradient since we
733 // frequently use it:
734 F_inv = invert(F);
735 tau = material->get_tau();
736 Jc = material->get_Jc();
737 dPsi_vol_dJ = material->get_dPsi_vol_dJ();
738 d2Psi_vol_dJ2 = material->get_d2Psi_vol_dJ2();
739 }
740
741 // We offer an interface to retrieve certain data. Here are the kinematic
742 // variables:
get_J_tilde() const743 double get_J_tilde() const
744 {
745 return material->get_J_tilde();
746 }
747
get_det_F() const748 double get_det_F() const
749 {
750 return material->get_det_F();
751 }
752
get_F_inv() const753 const Tensor<2, dim> &get_F_inv() const
754 {
755 return F_inv;
756 }
757
758 // ...and the kinetic variables. These are used in the material and
759 // global tangent matrix and residual assembly operations:
get_p_tilde() const760 double get_p_tilde() const
761 {
762 return material->get_p_tilde();
763 }
764
get_tau() const765 const SymmetricTensor<2, dim> &get_tau() const
766 {
767 return tau;
768 }
769
get_dPsi_vol_dJ() const770 double get_dPsi_vol_dJ() const
771 {
772 return dPsi_vol_dJ;
773 }
774
get_d2Psi_vol_dJ2() const775 double get_d2Psi_vol_dJ2() const
776 {
777 return d2Psi_vol_dJ2;
778 }
779
780 // And finally the tangent:
get_Jc() const781 const SymmetricTensor<4, dim> &get_Jc() const
782 {
783 return Jc;
784 }
785
786 // In terms of member functions, this class stores for the quadrature
787 // point it represents a copy of a material type in case different
788 // materials are used in different regions of the domain, as well as the
789 // inverse of the deformation gradient...
790 private:
791 std::shared_ptr<Material_Compressible_Neo_Hook_Three_Field<dim>> material;
792
793 Tensor<2, dim> F_inv;
794
795 // ... and stress-type variables along with the tangent $J\mathfrak{c}$:
796 SymmetricTensor<2, dim> tau;
797 double d2Psi_vol_dJ2;
798 double dPsi_vol_dJ;
799
800 SymmetricTensor<4, dim> Jc;
801 };
802
803
804 // @sect3{Quasi-static quasi-incompressible finite-strain solid}
805
806 // The Solid class is the central class in that it represents the problem at
807 // hand. It follows the usual scheme in that all it really has is a
808 // constructor, destructor and a <code>run()</code> function that dispatches
809 // all the work to private functions of this class:
810 template <int dim>
811 class Solid
812 {
813 public:
814 Solid(const std::string &input_file);
815
816 void run();
817
818 private:
819 // In the private section of this class, we first forward declare a number
820 // of objects that are used in parallelizing work using the WorkStream
821 // object (see the @ref threads module for more information on this).
822 //
823 // We declare such structures for the computation of tangent (stiffness)
824 // matrix, right hand side, static condensation, and for updating
825 // quadrature points:
826 struct PerTaskData_K;
827 struct ScratchData_K;
828
829 struct PerTaskData_RHS;
830 struct ScratchData_RHS;
831
832 struct PerTaskData_SC;
833 struct ScratchData_SC;
834
835 struct PerTaskData_UQPH;
836 struct ScratchData_UQPH;
837
838 // We start the collection of member functions with one that builds the
839 // grid:
840 void make_grid();
841
842 // Set up the finite element system to be solved:
843 void system_setup();
844
845 void determine_component_extractors();
846
847 // Several functions to assemble the system and right hand side matrices
848 // using multithreading. Each of them comes as a wrapper function, one
849 // that is executed to do the work in the WorkStream model on one cell,
850 // and one that copies the work done on this one cell into the global
851 // object that represents it:
852 void assemble_system_tangent();
853
854 void assemble_system_tangent_one_cell(
855 const typename DoFHandler<dim>::active_cell_iterator &cell,
856 ScratchData_K & scratch,
857 PerTaskData_K & data) const;
858
859 void copy_local_to_global_K(const PerTaskData_K &data);
860
861 void assemble_system_rhs();
862
863 void assemble_system_rhs_one_cell(
864 const typename DoFHandler<dim>::active_cell_iterator &cell,
865 ScratchData_RHS & scratch,
866 PerTaskData_RHS & data) const;
867
868 void copy_local_to_global_rhs(const PerTaskData_RHS &data);
869
870 void assemble_sc();
871
872 void assemble_sc_one_cell(
873 const typename DoFHandler<dim>::active_cell_iterator &cell,
874 ScratchData_SC & scratch,
875 PerTaskData_SC & data);
876
877 void copy_local_to_global_sc(const PerTaskData_SC &data);
878
879 // Apply Dirichlet boundary conditions on the displacement field
880 void make_constraints(const int &it_nr);
881
882 // Create and update the quadrature points. Here, no data needs to be
883 // copied into a global object, so the copy_local_to_global function is
884 // empty:
885 void setup_qph();
886
887 void update_qph_incremental(const BlockVector<double> &solution_delta);
888
889 void update_qph_incremental_one_cell(
890 const typename DoFHandler<dim>::active_cell_iterator &cell,
891 ScratchData_UQPH & scratch,
892 PerTaskData_UQPH & data);
893
copy_local_to_global_UQPH(const PerTaskData_UQPH &)894 void copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/)
895 {}
896
897 // Solve for the displacement using a Newton-Raphson method. We break this
898 // function into the nonlinear loop and the function that solves the
899 // linearized Newton-Raphson step:
900 void solve_nonlinear_timestep(BlockVector<double> &solution_delta);
901
902 std::pair<unsigned int, double>
903 solve_linear_system(BlockVector<double> &newton_update);
904
905 // Solution retrieval as well as post-processing and writing data to file:
906 BlockVector<double>
907 get_total_solution(const BlockVector<double> &solution_delta) const;
908
909 void output_results() const;
910
911 // Finally, some member variables that describe the current state: A
912 // collection of the parameters used to describe the problem setup...
913 Parameters::AllParameters parameters;
914
915 // ...the volume of the reference configuration...
916 double vol_reference;
917
918 // ...and description of the geometry on which the problem is solved:
919 Triangulation<dim> triangulation;
920
921 // Also, keep track of the current time and the time spent evaluating
922 // certain functions
923 Time time;
924 mutable TimerOutput timer;
925
926 // A storage object for quadrature point information. As opposed to
927 // step-18, deal.II's native quadrature point data manager is employed
928 // here.
929 CellDataStorage<typename Triangulation<dim>::cell_iterator,
930 PointHistory<dim>>
931 quadrature_point_history;
932
933 // A description of the finite-element system including the displacement
934 // polynomial degree, the degree-of-freedom handler, number of DoFs per
935 // cell and the extractor objects used to retrieve information from the
936 // solution vectors:
937 const unsigned int degree;
938 const FESystem<dim> fe;
939 DoFHandler<dim> dof_handler;
940 const unsigned int dofs_per_cell;
941 const FEValuesExtractors::Vector u_fe;
942 const FEValuesExtractors::Scalar p_fe;
943 const FEValuesExtractors::Scalar J_fe;
944
945 // Description of how the block-system is arranged. There are 3 blocks,
946 // the first contains a vector DOF $\mathbf{u}$ while the other two
947 // describe scalar DOFs, $\widetilde{p}$ and $\widetilde{J}$.
948 static const unsigned int n_blocks = 3;
949 static const unsigned int n_components = dim + 2;
950 static const unsigned int first_u_component = 0;
951 static const unsigned int p_component = dim;
952 static const unsigned int J_component = dim + 1;
953
954 enum
955 {
956 u_dof = 0,
957 p_dof = 1,
958 J_dof = 2
959 };
960
961 std::vector<types::global_dof_index> dofs_per_block;
962 std::vector<types::global_dof_index> element_indices_u;
963 std::vector<types::global_dof_index> element_indices_p;
964 std::vector<types::global_dof_index> element_indices_J;
965
966 // Rules for Gauss-quadrature on both the cell and faces. The number of
967 // quadrature points on both cells and faces is recorded.
968 const QGauss<dim> qf_cell;
969 const QGauss<dim - 1> qf_face;
970 const unsigned int n_q_points;
971 const unsigned int n_q_points_f;
972
973 // Objects that store the converged solution and right-hand side vectors,
974 // as well as the tangent matrix. There is an AffineConstraints object used
975 // to keep track of constraints. We make use of a sparsity pattern
976 // designed for a block system.
977 AffineConstraints<double> constraints;
978 BlockSparsityPattern sparsity_pattern;
979 BlockSparseMatrix<double> tangent_matrix;
980 BlockVector<double> system_rhs;
981 BlockVector<double> solution_n;
982
983 // Then define a number of variables to store norms and update norms and
984 // normalization factors.
985 struct Errors
986 {
ErrorsStep44::Solid::Errors987 Errors()
988 : norm(1.0)
989 , u(1.0)
990 , p(1.0)
991 , J(1.0)
992 {}
993
resetStep44::Solid::Errors994 void reset()
995 {
996 norm = 1.0;
997 u = 1.0;
998 p = 1.0;
999 J = 1.0;
1000 }
normalizeStep44::Solid::Errors1001 void normalize(const Errors &rhs)
1002 {
1003 if (rhs.norm != 0.0)
1004 norm /= rhs.norm;
1005 if (rhs.u != 0.0)
1006 u /= rhs.u;
1007 if (rhs.p != 0.0)
1008 p /= rhs.p;
1009 if (rhs.J != 0.0)
1010 J /= rhs.J;
1011 }
1012
1013 double norm, u, p, J;
1014 };
1015
1016 Errors error_residual, error_residual_0, error_residual_norm, error_update,
1017 error_update_0, error_update_norm;
1018
1019 // Methods to calculate error measures
1020 void get_error_residual(Errors &error_residual);
1021
1022 void get_error_update(const BlockVector<double> &newton_update,
1023 Errors & error_update);
1024
1025 std::pair<double, double> get_error_dilation() const;
1026
1027 // Compute the volume in the spatial configuration
1028 double compute_vol_current() const;
1029
1030 // Print information to screen in a pleasing way...
1031 static void print_conv_header();
1032
1033 void print_conv_footer();
1034 };
1035
1036 // @sect3{Implementation of the <code>Solid</code> class}
1037
1038 // @sect4{Public interface}
1039
1040 // We initialize the Solid class using data extracted from the parameter file.
1041 template <int dim>
Solid(const std::string & input_file)1042 Solid<dim>::Solid(const std::string &input_file)
1043 : parameters(input_file)
1044 , vol_reference(0.)
1045 , triangulation(Triangulation<dim>::maximum_smoothing)
1046 , time(parameters.end_time, parameters.delta_t)
1047 , timer(std::cout, TimerOutput::summary, TimerOutput::wall_times)
1048 , degree(parameters.poly_degree)
1049 ,
1050 // The Finite Element System is composed of dim continuous displacement
1051 // DOFs, and discontinuous pressure and dilatation DOFs. In an attempt to
1052 // satisfy the Babuska-Brezzi or LBB stability conditions (see Hughes
1053 // (2000)), we setup a $Q_n \times DGPM_{n-1} \times DGPM_{n-1}$
1054 // system. $Q_2 \times DGPM_1 \times DGPM_1$ elements satisfy this
1055 // condition, while $Q_1 \times DGPM_0 \times DGPM_0$ elements do
1056 // not. However, it has been shown that the latter demonstrate good
1057 // convergence characteristics nonetheless.
1058 fe(FE_Q<dim>(parameters.poly_degree),
1059 dim, // displacement
1060 FE_DGPMonomial<dim>(parameters.poly_degree - 1),
1061 1, // pressure
1062 FE_DGPMonomial<dim>(parameters.poly_degree - 1),
1063 1)
1064 , // dilatation
1065 dof_handler(triangulation)
1066 , dofs_per_cell(fe.n_dofs_per_cell())
1067 , u_fe(first_u_component)
1068 , p_fe(p_component)
1069 , J_fe(J_component)
1070 , dofs_per_block(n_blocks)
1071 , qf_cell(parameters.quad_order)
1072 , qf_face(parameters.quad_order)
1073 , n_q_points(qf_cell.size())
1074 , n_q_points_f(qf_face.size())
1075 {
1076 Assert(dim == 2 || dim == 3,
1077 ExcMessage("This problem only works in 2 or 3 space dimensions."));
1078 determine_component_extractors();
1079 }
1080
1081
1082 // In solving the quasi-static problem, the time becomes a loading parameter,
1083 // i.e. we increasing the loading linearly with time, making the two concepts
1084 // interchangeable. We choose to increment time linearly using a constant time
1085 // step size.
1086 //
1087 // We start the function with preprocessing, setting the initial dilatation
1088 // values, and then output the initial grid before starting the simulation
1089 // proper with the first time (and loading)
1090 // increment.
1091 //
1092 // Care must be taken (or at least some thought given) when imposing the
1093 // constraint $\widetilde{J}=1$ on the initial solution field. The constraint
1094 // corresponds to the determinant of the deformation gradient in the
1095 // undeformed configuration, which is the identity tensor. We use
1096 // FE_DGPMonomial bases to interpolate the dilatation field, thus we can't
1097 // simply set the corresponding dof to unity as they correspond to the
1098 // monomial coefficients. Thus we use the VectorTools::project function to do
1099 // the work for us. The VectorTools::project function requires an argument
1100 // indicating the hanging node constraints. We have none in this program
1101 // So we have to create a constraint object. In its original state, constraint
1102 // objects are unsorted, and have to be sorted (using the
1103 // AffineConstraints::close function) before they can be used. Have a look at
1104 // step-21 for more information. We only need to enforce the initial condition
1105 // on the dilatation. In order to do this, we make use of a
1106 // ComponentSelectFunction which acts as a mask and sets the J_component of
1107 // n_components to 1. This is exactly what we want. Have a look at its usage
1108 // in step-20 for more information.
1109 template <int dim>
run()1110 void Solid<dim>::run()
1111 {
1112 make_grid();
1113 system_setup();
1114 {
1115 AffineConstraints<double> constraints;
1116 constraints.close();
1117
1118 const ComponentSelectFunction<dim> J_mask(J_component, n_components);
1119
1120 VectorTools::project(
1121 dof_handler, constraints, QGauss<dim>(degree + 2), J_mask, solution_n);
1122 }
1123 output_results();
1124 time.increment();
1125
1126 // We then declare the incremental solution update $\varDelta
1127 // \mathbf{\Xi} \dealcoloneq \{\varDelta \mathbf{u},\varDelta \widetilde{p},
1128 // \varDelta \widetilde{J} \}$ and start the loop over the time domain.
1129 //
1130 // At the beginning, we reset the solution update for this time step...
1131 BlockVector<double> solution_delta(dofs_per_block);
1132 while (time.current() < time.end())
1133 {
1134 solution_delta = 0.0;
1135
1136 // ...solve the current time step and update total solution vector
1137 // $\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} +
1138 // \varDelta \mathbf{\Xi}$...
1139 solve_nonlinear_timestep(solution_delta);
1140 solution_n += solution_delta;
1141
1142 // ...and plot the results before moving on happily to the next time
1143 // step:
1144 output_results();
1145 time.increment();
1146 }
1147 }
1148
1149
1150 // @sect3{Private interface}
1151
1152 // @sect4{Threading-building-blocks structures}
1153
1154 // The first group of private member functions is related to parallelization.
1155 // We use the Threading Building Blocks library (TBB) to perform as many
1156 // computationally intensive distributed tasks as possible. In particular, we
1157 // assemble the tangent matrix and right hand side vector, the static
1158 // condensation contributions, and update data stored at the quadrature points
1159 // using TBB. Our main tool for this is the WorkStream class (see the @ref
1160 // threads module for more information).
1161
1162 // Firstly we deal with the tangent matrix assembly structures. The
1163 // PerTaskData object stores local contributions.
1164 template <int dim>
1165 struct Solid<dim>::PerTaskData_K
1166 {
1167 FullMatrix<double> cell_matrix;
1168 std::vector<types::global_dof_index> local_dof_indices;
1169
PerTaskData_KStep44::Solid::PerTaskData_K1170 PerTaskData_K(const unsigned int dofs_per_cell)
1171 : cell_matrix(dofs_per_cell, dofs_per_cell)
1172 , local_dof_indices(dofs_per_cell)
1173 {}
1174
resetStep44::Solid::PerTaskData_K1175 void reset()
1176 {
1177 cell_matrix = 0.0;
1178 }
1179 };
1180
1181
1182 // On the other hand, the ScratchData object stores the larger objects such as
1183 // the shape-function values array (<code>Nx</code>) and a shape function
1184 // gradient and symmetric gradient vector which we will use during the
1185 // assembly.
1186 template <int dim>
1187 struct Solid<dim>::ScratchData_K
1188 {
1189 FEValues<dim> fe_values;
1190
1191 std::vector<std::vector<double>> Nx;
1192 std::vector<std::vector<Tensor<2, dim>>> grad_Nx;
1193 std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
1194
ScratchData_KStep44::Solid::ScratchData_K1195 ScratchData_K(const FiniteElement<dim> &fe_cell,
1196 const QGauss<dim> & qf_cell,
1197 const UpdateFlags uf_cell)
1198 : fe_values(fe_cell, qf_cell, uf_cell)
1199 , Nx(qf_cell.size(), std::vector<double>(fe_cell.n_dofs_per_cell()))
1200 , grad_Nx(qf_cell.size(),
1201 std::vector<Tensor<2, dim>>(fe_cell.n_dofs_per_cell()))
1202 , symm_grad_Nx(qf_cell.size(),
1203 std::vector<SymmetricTensor<2, dim>>(
1204 fe_cell.n_dofs_per_cell()))
1205 {}
1206
ScratchData_KStep44::Solid::ScratchData_K1207 ScratchData_K(const ScratchData_K &rhs)
1208 : fe_values(rhs.fe_values.get_fe(),
1209 rhs.fe_values.get_quadrature(),
1210 rhs.fe_values.get_update_flags())
1211 , Nx(rhs.Nx)
1212 , grad_Nx(rhs.grad_Nx)
1213 , symm_grad_Nx(rhs.symm_grad_Nx)
1214 {}
1215
resetStep44::Solid::ScratchData_K1216 void reset()
1217 {
1218 const unsigned int n_q_points = Nx.size();
1219 const unsigned int n_dofs_per_cell = Nx[0].size();
1220 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1221 {
1222 Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
1223 Assert(grad_Nx[q_point].size() == n_dofs_per_cell,
1224 ExcInternalError());
1225 Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
1226 ExcInternalError());
1227 for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
1228 {
1229 Nx[q_point][k] = 0.0;
1230 grad_Nx[q_point][k] = 0.0;
1231 symm_grad_Nx[q_point][k] = 0.0;
1232 }
1233 }
1234 }
1235 };
1236
1237 // Next, the same approach is used for the right-hand side assembly. The
1238 // PerTaskData object again stores local contributions and the ScratchData
1239 // object the shape function object and precomputed values vector:
1240 template <int dim>
1241 struct Solid<dim>::PerTaskData_RHS
1242 {
1243 Vector<double> cell_rhs;
1244 std::vector<types::global_dof_index> local_dof_indices;
1245
PerTaskData_RHSStep44::Solid::PerTaskData_RHS1246 PerTaskData_RHS(const unsigned int dofs_per_cell)
1247 : cell_rhs(dofs_per_cell)
1248 , local_dof_indices(dofs_per_cell)
1249 {}
1250
resetStep44::Solid::PerTaskData_RHS1251 void reset()
1252 {
1253 cell_rhs = 0.0;
1254 }
1255 };
1256
1257
1258 template <int dim>
1259 struct Solid<dim>::ScratchData_RHS
1260 {
1261 FEValues<dim> fe_values;
1262 FEFaceValues<dim> fe_face_values;
1263
1264 std::vector<std::vector<double>> Nx;
1265 std::vector<std::vector<SymmetricTensor<2, dim>>> symm_grad_Nx;
1266
ScratchData_RHSStep44::Solid::ScratchData_RHS1267 ScratchData_RHS(const FiniteElement<dim> &fe_cell,
1268 const QGauss<dim> & qf_cell,
1269 const UpdateFlags uf_cell,
1270 const QGauss<dim - 1> & qf_face,
1271 const UpdateFlags uf_face)
1272 : fe_values(fe_cell, qf_cell, uf_cell)
1273 , fe_face_values(fe_cell, qf_face, uf_face)
1274 , Nx(qf_cell.size(), std::vector<double>(fe_cell.n_dofs_per_cell()))
1275 , symm_grad_Nx(qf_cell.size(),
1276 std::vector<SymmetricTensor<2, dim>>(
1277 fe_cell.n_dofs_per_cell()))
1278 {}
1279
ScratchData_RHSStep44::Solid::ScratchData_RHS1280 ScratchData_RHS(const ScratchData_RHS &rhs)
1281 : fe_values(rhs.fe_values.get_fe(),
1282 rhs.fe_values.get_quadrature(),
1283 rhs.fe_values.get_update_flags())
1284 , fe_face_values(rhs.fe_face_values.get_fe(),
1285 rhs.fe_face_values.get_quadrature(),
1286 rhs.fe_face_values.get_update_flags())
1287 , Nx(rhs.Nx)
1288 , symm_grad_Nx(rhs.symm_grad_Nx)
1289 {}
1290
resetStep44::Solid::ScratchData_RHS1291 void reset()
1292 {
1293 const unsigned int n_q_points = Nx.size();
1294 const unsigned int n_dofs_per_cell = Nx[0].size();
1295 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1296 {
1297 Assert(Nx[q_point].size() == n_dofs_per_cell, ExcInternalError());
1298 Assert(symm_grad_Nx[q_point].size() == n_dofs_per_cell,
1299 ExcInternalError());
1300 for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
1301 {
1302 Nx[q_point][k] = 0.0;
1303 symm_grad_Nx[q_point][k] = 0.0;
1304 }
1305 }
1306 }
1307 };
1308
1309 // Then we define structures to assemble the statically condensed tangent
1310 // matrix. Recall that we wish to solve for a displacement-based formulation.
1311 // We do the condensation at the element level as the $\widetilde{p}$ and
1312 // $\widetilde{J}$ fields are element-wise discontinuous. As these operations
1313 // are matrix-based, we need to setup a number of matrices to store the local
1314 // contributions from a number of the tangent matrix sub-blocks. We place
1315 // these in the PerTaskData struct.
1316 //
1317 // We choose not to reset any data in the <code>reset()</code> function as the
1318 // matrix extraction and replacement tools will take care of this
1319 template <int dim>
1320 struct Solid<dim>::PerTaskData_SC
1321 {
1322 FullMatrix<double> cell_matrix;
1323 std::vector<types::global_dof_index> local_dof_indices;
1324
1325 FullMatrix<double> k_orig;
1326 FullMatrix<double> k_pu;
1327 FullMatrix<double> k_pJ;
1328 FullMatrix<double> k_JJ;
1329 FullMatrix<double> k_pJ_inv;
1330 FullMatrix<double> k_bbar;
1331 FullMatrix<double> A;
1332 FullMatrix<double> B;
1333 FullMatrix<double> C;
1334
PerTaskData_SCStep44::Solid::PerTaskData_SC1335 PerTaskData_SC(const unsigned int dofs_per_cell,
1336 const unsigned int n_u,
1337 const unsigned int n_p,
1338 const unsigned int n_J)
1339 : cell_matrix(dofs_per_cell, dofs_per_cell)
1340 , local_dof_indices(dofs_per_cell)
1341 , k_orig(dofs_per_cell, dofs_per_cell)
1342 , k_pu(n_p, n_u)
1343 , k_pJ(n_p, n_J)
1344 , k_JJ(n_J, n_J)
1345 , k_pJ_inv(n_p, n_J)
1346 , k_bbar(n_u, n_u)
1347 , A(n_J, n_u)
1348 , B(n_J, n_u)
1349 , C(n_p, n_u)
1350 {}
1351
resetStep44::Solid::PerTaskData_SC1352 void reset()
1353 {}
1354 };
1355
1356
1357 // The ScratchData object for the operations we wish to perform here is empty
1358 // since we need no temporary data, but it still needs to be defined for the
1359 // current implementation of TBB in deal.II. So we create a dummy struct for
1360 // this purpose.
1361 template <int dim>
1362 struct Solid<dim>::ScratchData_SC
1363 {
resetStep44::Solid::ScratchData_SC1364 void reset()
1365 {}
1366 };
1367
1368
1369 // And finally we define the structures to assist with updating the quadrature
1370 // point information. Similar to the SC assembly process, we do not need the
1371 // PerTaskData object (since there is nothing to store here) but must define
1372 // one nonetheless. Note that this is because for the operation that we have
1373 // here -- updating the data on quadrature points -- the operation is purely
1374 // local: the things we do on every cell get consumed on every cell, without
1375 // any global aggregation operation as is usually the case when using the
1376 // WorkStream class. The fact that we still have to define a per-task data
1377 // structure points to the fact that the WorkStream class may be ill-suited to
1378 // this operation (we could, in principle simply create a new task using
1379 // Threads::new_task for each cell) but there is not much harm done to doing
1380 // it this way anyway.
1381 // Furthermore, should there be different material models associated with a
1382 // quadrature point, requiring varying levels of computational expense, then
1383 // the method used here could be advantageous.
1384 template <int dim>
1385 struct Solid<dim>::PerTaskData_UQPH
1386 {
resetStep44::Solid::PerTaskData_UQPH1387 void reset()
1388 {}
1389 };
1390
1391
1392 // The ScratchData object will be used to store an alias for the solution
1393 // vector so that we don't have to copy this large data structure. We then
1394 // define a number of vectors to extract the solution values and gradients at
1395 // the quadrature points.
1396 template <int dim>
1397 struct Solid<dim>::ScratchData_UQPH
1398 {
1399 const BlockVector<double> &solution_total;
1400
1401 std::vector<Tensor<2, dim>> solution_grads_u_total;
1402 std::vector<double> solution_values_p_total;
1403 std::vector<double> solution_values_J_total;
1404
1405 FEValues<dim> fe_values;
1406
ScratchData_UQPHStep44::Solid::ScratchData_UQPH1407 ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
1408 const QGauss<dim> & qf_cell,
1409 const UpdateFlags uf_cell,
1410 const BlockVector<double> &solution_total)
1411 : solution_total(solution_total)
1412 , solution_grads_u_total(qf_cell.size())
1413 , solution_values_p_total(qf_cell.size())
1414 , solution_values_J_total(qf_cell.size())
1415 , fe_values(fe_cell, qf_cell, uf_cell)
1416 {}
1417
ScratchData_UQPHStep44::Solid::ScratchData_UQPH1418 ScratchData_UQPH(const ScratchData_UQPH &rhs)
1419 : solution_total(rhs.solution_total)
1420 , solution_grads_u_total(rhs.solution_grads_u_total)
1421 , solution_values_p_total(rhs.solution_values_p_total)
1422 , solution_values_J_total(rhs.solution_values_J_total)
1423 , fe_values(rhs.fe_values.get_fe(),
1424 rhs.fe_values.get_quadrature(),
1425 rhs.fe_values.get_update_flags())
1426 {}
1427
resetStep44::Solid::ScratchData_UQPH1428 void reset()
1429 {
1430 const unsigned int n_q_points = solution_grads_u_total.size();
1431 for (unsigned int q = 0; q < n_q_points; ++q)
1432 {
1433 solution_grads_u_total[q] = 0.0;
1434 solution_values_p_total[q] = 0.0;
1435 solution_values_J_total[q] = 0.0;
1436 }
1437 }
1438 };
1439
1440
1441 // @sect4{Solid::make_grid}
1442
1443 // On to the first of the private member functions. Here we create the
1444 // triangulation of the domain, for which we choose the scaled cube with each
1445 // face given a boundary ID number. The grid must be refined at least once
1446 // for the indentation problem.
1447 //
1448 // We then determine the volume of the reference configuration and print it
1449 // for comparison:
1450 template <int dim>
make_grid()1451 void Solid<dim>::make_grid()
1452 {
1453 GridGenerator::hyper_rectangle(
1454 triangulation,
1455 (dim == 3 ? Point<dim>(0.0, 0.0, 0.0) : Point<dim>(0.0, 0.0)),
1456 (dim == 3 ? Point<dim>(1.0, 1.0, 1.0) : Point<dim>(1.0, 1.0)),
1457 true);
1458 GridTools::scale(parameters.scale, triangulation);
1459 triangulation.refine_global(std::max(1U, parameters.global_refinement));
1460
1461 vol_reference = GridTools::volume(triangulation);
1462 std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl;
1463
1464 // Since we wish to apply a Neumann BC to a patch on the top surface, we
1465 // must find the cell faces in this part of the domain and mark them with
1466 // a distinct boundary ID number. The faces we are looking for are on the
1467 // +y surface and will get boundary ID 6 (zero through five are already
1468 // used when creating the six faces of the cube domain):
1469 for (const auto &cell : triangulation.active_cell_iterators())
1470 for (const auto &face : cell->face_iterators())
1471 {
1472 if (face->at_boundary() == true &&
1473 face->center()[1] == 1.0 * parameters.scale)
1474 {
1475 if (dim == 3)
1476 {
1477 if (face->center()[0] < 0.5 * parameters.scale &&
1478 face->center()[2] < 0.5 * parameters.scale)
1479 face->set_boundary_id(6);
1480 }
1481 else
1482 {
1483 if (face->center()[0] < 0.5 * parameters.scale)
1484 face->set_boundary_id(6);
1485 }
1486 }
1487 }
1488 }
1489
1490
1491 // @sect4{Solid::system_setup}
1492
1493 // Next we describe how the FE system is setup. We first determine the number
1494 // of components per block. Since the displacement is a vector component, the
1495 // first dim components belong to it, while the next two describe scalar
1496 // pressure and dilatation DOFs.
1497 template <int dim>
system_setup()1498 void Solid<dim>::system_setup()
1499 {
1500 timer.enter_subsection("Setup system");
1501
1502 std::vector<unsigned int> block_component(n_components,
1503 u_dof); // Displacement
1504 block_component[p_component] = p_dof; // Pressure
1505 block_component[J_component] = J_dof; // Dilatation
1506
1507 // The DOF handler is then initialized and we renumber the grid in an
1508 // efficient manner. We also record the number of DOFs per block.
1509 dof_handler.distribute_dofs(fe);
1510 DoFRenumbering::Cuthill_McKee(dof_handler);
1511 DoFRenumbering::component_wise(dof_handler, block_component);
1512
1513 dofs_per_block =
1514 DoFTools::count_dofs_per_fe_block(dof_handler, block_component);
1515
1516 std::cout << "Triangulation:"
1517 << "\n\t Number of active cells: "
1518 << triangulation.n_active_cells()
1519 << "\n\t Number of degrees of freedom: " << dof_handler.n_dofs()
1520 << std::endl;
1521
1522 // Setup the sparsity pattern and tangent matrix
1523 tangent_matrix.clear();
1524 {
1525 const types::global_dof_index n_dofs_u = dofs_per_block[u_dof];
1526 const types::global_dof_index n_dofs_p = dofs_per_block[p_dof];
1527 const types::global_dof_index n_dofs_J = dofs_per_block[J_dof];
1528
1529 BlockDynamicSparsityPattern dsp(n_blocks, n_blocks);
1530
1531 dsp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u);
1532 dsp.block(u_dof, p_dof).reinit(n_dofs_u, n_dofs_p);
1533 dsp.block(u_dof, J_dof).reinit(n_dofs_u, n_dofs_J);
1534
1535 dsp.block(p_dof, u_dof).reinit(n_dofs_p, n_dofs_u);
1536 dsp.block(p_dof, p_dof).reinit(n_dofs_p, n_dofs_p);
1537 dsp.block(p_dof, J_dof).reinit(n_dofs_p, n_dofs_J);
1538
1539 dsp.block(J_dof, u_dof).reinit(n_dofs_J, n_dofs_u);
1540 dsp.block(J_dof, p_dof).reinit(n_dofs_J, n_dofs_p);
1541 dsp.block(J_dof, J_dof).reinit(n_dofs_J, n_dofs_J);
1542 dsp.collect_sizes();
1543
1544 // The global system matrix initially has the following structure
1545 // @f{align*}
1546 // \underbrace{\begin{bmatrix}
1547 // \mathsf{\mathbf{K}}_{uu} & \mathsf{\mathbf{K}}_{u\widetilde{p}} &
1548 // \mathbf{0}
1549 // \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
1550 // \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}
1551 // \\ \mathbf{0} & \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
1552 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
1553 // \end{bmatrix}}_{\mathsf{\mathbf{K}}(\mathbf{\Xi}_{\textrm{i}})}
1554 // \underbrace{\begin{bmatrix}
1555 // d \mathsf{u}
1556 // \\ d \widetilde{\mathsf{\mathbf{p}}}
1557 // \\ d \widetilde{\mathsf{\mathbf{J}}}
1558 // \end{bmatrix}}_{d \mathbf{\Xi}}
1559 // =
1560 // \underbrace{\begin{bmatrix}
1561 // \mathsf{\mathbf{F}}_{u}(\mathbf{u}_{\textrm{i}})
1562 // \\ \mathsf{\mathbf{F}}_{\widetilde{p}}(\widetilde{p}_{\textrm{i}})
1563 // \\ \mathsf{\mathbf{F}}_{\widetilde{J}}(\widetilde{J}_{\textrm{i}})
1564 //\end{bmatrix}}_{ \mathsf{\mathbf{F}}(\mathbf{\Xi}_{\textrm{i}}) } \, .
1565 // @f}
1566 // We optimize the sparsity pattern to reflect this structure
1567 // and prevent unnecessary data creation for the right-diagonal
1568 // block components.
1569 Table<2, DoFTools::Coupling> coupling(n_components, n_components);
1570 for (unsigned int ii = 0; ii < n_components; ++ii)
1571 for (unsigned int jj = 0; jj < n_components; ++jj)
1572 if (((ii < p_component) && (jj == J_component)) ||
1573 ((ii == J_component) && (jj < p_component)) ||
1574 ((ii == p_component) && (jj == p_component)))
1575 coupling[ii][jj] = DoFTools::none;
1576 else
1577 coupling[ii][jj] = DoFTools::always;
1578 DoFTools::make_sparsity_pattern(
1579 dof_handler, coupling, dsp, constraints, false);
1580 sparsity_pattern.copy_from(dsp);
1581 }
1582
1583 tangent_matrix.reinit(sparsity_pattern);
1584
1585 // We then set up storage vectors
1586 system_rhs.reinit(dofs_per_block);
1587 system_rhs.collect_sizes();
1588
1589 solution_n.reinit(dofs_per_block);
1590 solution_n.collect_sizes();
1591
1592 // ...and finally set up the quadrature
1593 // point history:
1594 setup_qph();
1595
1596 timer.leave_subsection();
1597 }
1598
1599
1600 // @sect4{Solid::determine_component_extractors}
1601 // Next we compute some information from the FE system that describes which
1602 // local element DOFs are attached to which block component. This is used
1603 // later to extract sub-blocks from the global matrix.
1604 //
1605 // In essence, all we need is for the FESystem object to indicate to which
1606 // block component a DOF on the reference cell is attached to. Currently, the
1607 // interpolation fields are setup such that 0 indicates a displacement DOF, 1
1608 // a pressure DOF and 2 a dilatation DOF.
1609 template <int dim>
determine_component_extractors()1610 void Solid<dim>::determine_component_extractors()
1611 {
1612 element_indices_u.clear();
1613 element_indices_p.clear();
1614 element_indices_J.clear();
1615
1616 for (unsigned int k = 0; k < fe.n_dofs_per_cell(); ++k)
1617 {
1618 const unsigned int k_group = fe.system_to_base_index(k).first.first;
1619 if (k_group == u_dof)
1620 element_indices_u.push_back(k);
1621 else if (k_group == p_dof)
1622 element_indices_p.push_back(k);
1623 else if (k_group == J_dof)
1624 element_indices_J.push_back(k);
1625 else
1626 {
1627 Assert(k_group <= J_dof, ExcInternalError());
1628 }
1629 }
1630 }
1631
1632 // @sect4{Solid::setup_qph}
1633 // The method used to store quadrature information is already described in
1634 // step-18. Here we implement a similar setup for a SMP machine.
1635 //
1636 // Firstly the actual QPH data objects are created. This must be done only
1637 // once the grid is refined to its finest level.
1638 template <int dim>
setup_qph()1639 void Solid<dim>::setup_qph()
1640 {
1641 std::cout << " Setting up quadrature point data..." << std::endl;
1642
1643 quadrature_point_history.initialize(triangulation.begin_active(),
1644 triangulation.end(),
1645 n_q_points);
1646
1647 // Next we setup the initial quadrature point data.
1648 // Note that when the quadrature point data is retrieved,
1649 // it is returned as a vector of smart pointers.
1650 for (const auto &cell : triangulation.active_cell_iterators())
1651 {
1652 const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
1653 quadrature_point_history.get_data(cell);
1654 Assert(lqph.size() == n_q_points, ExcInternalError());
1655
1656 for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
1657 lqph[q_point]->setup_lqp(parameters);
1658 }
1659 }
1660
1661 // @sect4{Solid::update_qph_incremental}
1662 // As the update of QP information occurs frequently and involves a number of
1663 // expensive operations, we define a multithreaded approach to distributing
1664 // the task across a number of CPU cores.
1665 //
1666 // To start this, we first we need to obtain the total solution as it stands
1667 // at this Newton increment and then create the initial copy of the scratch
1668 // and copy data objects:
1669 template <int dim>
1670 void
update_qph_incremental(const BlockVector<double> & solution_delta)1671 Solid<dim>::update_qph_incremental(const BlockVector<double> &solution_delta)
1672 {
1673 timer.enter_subsection("Update QPH data");
1674 std::cout << " UQPH " << std::flush;
1675
1676 const BlockVector<double> solution_total(
1677 get_total_solution(solution_delta));
1678
1679 const UpdateFlags uf_UQPH(update_values | update_gradients);
1680 PerTaskData_UQPH per_task_data_UQPH;
1681 ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total);
1682
1683 // We then pass them and the one-cell update function to the WorkStream to
1684 // be processed:
1685 WorkStream::run(dof_handler.active_cell_iterators(),
1686 *this,
1687 &Solid::update_qph_incremental_one_cell,
1688 &Solid::copy_local_to_global_UQPH,
1689 scratch_data_UQPH,
1690 per_task_data_UQPH);
1691
1692 timer.leave_subsection();
1693 }
1694
1695
1696 // Now we describe how we extract data from the solution vector and pass it
1697 // along to each QP storage object for processing.
1698 template <int dim>
update_qph_incremental_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,ScratchData_UQPH & scratch,PerTaskData_UQPH &)1699 void Solid<dim>::update_qph_incremental_one_cell(
1700 const typename DoFHandler<dim>::active_cell_iterator &cell,
1701 ScratchData_UQPH & scratch,
1702 PerTaskData_UQPH & /*data*/)
1703 {
1704 const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
1705 quadrature_point_history.get_data(cell);
1706 Assert(lqph.size() == n_q_points, ExcInternalError());
1707
1708 Assert(scratch.solution_grads_u_total.size() == n_q_points,
1709 ExcInternalError());
1710 Assert(scratch.solution_values_p_total.size() == n_q_points,
1711 ExcInternalError());
1712 Assert(scratch.solution_values_J_total.size() == n_q_points,
1713 ExcInternalError());
1714
1715 scratch.reset();
1716
1717 // We first need to find the values and gradients at quadrature points
1718 // inside the current cell and then we update each local QP using the
1719 // displacement gradient and total pressure and dilatation solution
1720 // values:
1721 scratch.fe_values.reinit(cell);
1722 scratch.fe_values[u_fe].get_function_gradients(
1723 scratch.solution_total, scratch.solution_grads_u_total);
1724 scratch.fe_values[p_fe].get_function_values(
1725 scratch.solution_total, scratch.solution_values_p_total);
1726 scratch.fe_values[J_fe].get_function_values(
1727 scratch.solution_total, scratch.solution_values_J_total);
1728
1729 for (const unsigned int q_point :
1730 scratch.fe_values.quadrature_point_indices())
1731 lqph[q_point]->update_values(scratch.solution_grads_u_total[q_point],
1732 scratch.solution_values_p_total[q_point],
1733 scratch.solution_values_J_total[q_point]);
1734 }
1735
1736
1737 // @sect4{Solid::solve_nonlinear_timestep}
1738
1739 // The next function is the driver method for the Newton-Raphson scheme. At
1740 // its top we create a new vector to store the current Newton update step,
1741 // reset the error storage objects and print solver header.
1742 template <int dim>
solve_nonlinear_timestep(BlockVector<double> & solution_delta)1743 void Solid<dim>::solve_nonlinear_timestep(BlockVector<double> &solution_delta)
1744 {
1745 std::cout << std::endl
1746 << "Timestep " << time.get_timestep() << " @ " << time.current()
1747 << "s" << std::endl;
1748
1749 BlockVector<double> newton_update(dofs_per_block);
1750
1751 error_residual.reset();
1752 error_residual_0.reset();
1753 error_residual_norm.reset();
1754 error_update.reset();
1755 error_update_0.reset();
1756 error_update_norm.reset();
1757
1758 print_conv_header();
1759
1760 // We now perform a number of Newton iterations to iteratively solve the
1761 // nonlinear problem. Since the problem is fully nonlinear and we are
1762 // using a full Newton method, the data stored in the tangent matrix and
1763 // right-hand side vector is not reusable and must be cleared at each
1764 // Newton step. We then initially build the right-hand side vector to
1765 // check for convergence (and store this value in the first iteration).
1766 // The unconstrained DOFs of the rhs vector hold the out-of-balance
1767 // forces. The building is done before assembling the system matrix as the
1768 // latter is an expensive operation and we can potentially avoid an extra
1769 // assembly process by not assembling the tangent matrix when convergence
1770 // is attained.
1771 unsigned int newton_iteration = 0;
1772 for (; newton_iteration < parameters.max_iterations_NR; ++newton_iteration)
1773 {
1774 std::cout << " " << std::setw(2) << newton_iteration << " "
1775 << std::flush;
1776
1777 tangent_matrix = 0.0;
1778 system_rhs = 0.0;
1779
1780 assemble_system_rhs();
1781 get_error_residual(error_residual);
1782
1783 if (newton_iteration == 0)
1784 error_residual_0 = error_residual;
1785
1786 // We can now determine the normalized residual error and check for
1787 // solution convergence:
1788 error_residual_norm = error_residual;
1789 error_residual_norm.normalize(error_residual_0);
1790
1791 if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u &&
1792 error_residual_norm.u <= parameters.tol_f)
1793 {
1794 std::cout << " CONVERGED! " << std::endl;
1795 print_conv_footer();
1796
1797 break;
1798 }
1799
1800 // If we have decided that we want to continue with the iteration, we
1801 // assemble the tangent, make and impose the Dirichlet constraints,
1802 // and do the solve of the linearised system:
1803 assemble_system_tangent();
1804 make_constraints(newton_iteration);
1805 constraints.condense(tangent_matrix, system_rhs);
1806
1807 const std::pair<unsigned int, double> lin_solver_output =
1808 solve_linear_system(newton_update);
1809
1810 get_error_update(newton_update, error_update);
1811 if (newton_iteration == 0)
1812 error_update_0 = error_update;
1813
1814 // We can now determine the normalized Newton update error, and
1815 // perform the actual update of the solution increment for the current
1816 // time step, update all quadrature point information pertaining to
1817 // this new displacement and stress state and continue iterating:
1818 error_update_norm = error_update;
1819 error_update_norm.normalize(error_update_0);
1820
1821 solution_delta += newton_update;
1822 update_qph_incremental(solution_delta);
1823
1824 std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7)
1825 << std::scientific << lin_solver_output.first << " "
1826 << lin_solver_output.second << " "
1827 << error_residual_norm.norm << " " << error_residual_norm.u
1828 << " " << error_residual_norm.p << " "
1829 << error_residual_norm.J << " " << error_update_norm.norm
1830 << " " << error_update_norm.u << " " << error_update_norm.p
1831 << " " << error_update_norm.J << " " << std::endl;
1832 }
1833
1834 // At the end, if it turns out that we have in fact done more iterations
1835 // than the parameter file allowed, we raise an exception that can be
1836 // caught in the main() function. The call <code>AssertThrow(condition,
1837 // exc_object)</code> is in essence equivalent to <code>if (!cond) throw
1838 // exc_object;</code> but the former form fills certain fields in the
1839 // exception object that identify the location (filename and line number)
1840 // where the exception was raised to make it simpler to identify where the
1841 // problem happened.
1842 AssertThrow(newton_iteration < parameters.max_iterations_NR,
1843 ExcMessage("No convergence in nonlinear solver!"));
1844 }
1845
1846
1847 // @sect4{Solid::print_conv_header and Solid::print_conv_footer}
1848
1849 // This program prints out data in a nice table that is updated
1850 // on a per-iteration basis. The next two functions set up the table
1851 // header and footer:
1852 template <int dim>
print_conv_header()1853 void Solid<dim>::print_conv_header()
1854 {
1855 static const unsigned int l_width = 155;
1856
1857 for (unsigned int i = 0; i < l_width; ++i)
1858 std::cout << "_";
1859 std::cout << std::endl;
1860
1861 std::cout << " SOLVER STEP "
1862 << " | LIN_IT LIN_RES RES_NORM "
1863 << " RES_U RES_P RES_J NU_NORM "
1864 << " NU_U NU_P NU_J " << std::endl;
1865
1866 for (unsigned int i = 0; i < l_width; ++i)
1867 std::cout << "_";
1868 std::cout << std::endl;
1869 }
1870
1871
1872
1873 template <int dim>
print_conv_footer()1874 void Solid<dim>::print_conv_footer()
1875 {
1876 static const unsigned int l_width = 155;
1877
1878 for (unsigned int i = 0; i < l_width; ++i)
1879 std::cout << "_";
1880 std::cout << std::endl;
1881
1882 const std::pair<double, double> error_dil = get_error_dilation();
1883
1884 std::cout << "Relative errors:" << std::endl
1885 << "Displacement:\t" << error_update.u / error_update_0.u
1886 << std::endl
1887 << "Force: \t\t" << error_residual.u / error_residual_0.u
1888 << std::endl
1889 << "Dilatation:\t" << error_dil.first << std::endl
1890 << "v / V_0:\t" << error_dil.second * vol_reference << " / "
1891 << vol_reference << " = " << error_dil.second << std::endl;
1892 }
1893
1894
1895 // @sect4{Solid::get_error_dilation}
1896
1897 // Calculate the volume of the domain in the spatial configuration
1898 template <int dim>
compute_vol_current() const1899 double Solid<dim>::compute_vol_current() const
1900 {
1901 double vol_current = 0.0;
1902
1903 FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
1904
1905 for (const auto &cell : triangulation.active_cell_iterators())
1906 {
1907 fe_values.reinit(cell);
1908
1909 // In contrast to that which was previously called for,
1910 // in this instance the quadrature point data is specifically
1911 // non-modifiable since we will only be accessing data.
1912 // We ensure that the right get_data function is called by
1913 // marking this update function as constant.
1914 const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
1915 quadrature_point_history.get_data(cell);
1916 Assert(lqph.size() == n_q_points, ExcInternalError());
1917
1918 for (const unsigned int q_point : fe_values.quadrature_point_indices())
1919 {
1920 const double det_F_qp = lqph[q_point]->get_det_F();
1921 const double JxW = fe_values.JxW(q_point);
1922
1923 vol_current += det_F_qp * JxW;
1924 }
1925 }
1926 Assert(vol_current > 0.0, ExcInternalError());
1927 return vol_current;
1928 }
1929
1930 // Calculate how well the dilatation $\widetilde{J}$ agrees with $J
1931 // \dealcoloneq \textrm{det}\ \mathbf{F}$ from the $L^2$ error $ \bigl[
1932 // \int_{\Omega_0} {[ J - \widetilde{J}]}^{2}\textrm{d}V \bigr]^{1/2}$.
1933 // We also return the ratio of the current volume of the
1934 // domain to the reference volume. This is of interest for incompressible
1935 // media where we want to check how well the isochoric constraint has been
1936 // enforced.
1937 template <int dim>
get_error_dilation() const1938 std::pair<double, double> Solid<dim>::get_error_dilation() const
1939 {
1940 double dil_L2_error = 0.0;
1941
1942 FEValues<dim> fe_values(fe, qf_cell, update_JxW_values);
1943
1944 for (const auto &cell : triangulation.active_cell_iterators())
1945 {
1946 fe_values.reinit(cell);
1947
1948 const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
1949 quadrature_point_history.get_data(cell);
1950 Assert(lqph.size() == n_q_points, ExcInternalError());
1951
1952 for (const unsigned int q_point : fe_values.quadrature_point_indices())
1953 {
1954 const double det_F_qp = lqph[q_point]->get_det_F();
1955 const double J_tilde_qp = lqph[q_point]->get_J_tilde();
1956 const double the_error_qp_squared =
1957 std::pow((det_F_qp - J_tilde_qp), 2);
1958 const double JxW = fe_values.JxW(q_point);
1959
1960 dil_L2_error += the_error_qp_squared * JxW;
1961 }
1962 }
1963
1964 return std::make_pair(std::sqrt(dil_L2_error),
1965 compute_vol_current() / vol_reference);
1966 }
1967
1968
1969 // @sect4{Solid::get_error_residual}
1970
1971 // Determine the true residual error for the problem. That is, determine the
1972 // error in the residual for the unconstrained degrees of freedom. Note that
1973 // to do so, we need to ignore constrained DOFs by setting the residual in
1974 // these vector components to zero.
1975 template <int dim>
get_error_residual(Errors & error_residual)1976 void Solid<dim>::get_error_residual(Errors &error_residual)
1977 {
1978 BlockVector<double> error_res(dofs_per_block);
1979
1980 for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
1981 if (!constraints.is_constrained(i))
1982 error_res(i) = system_rhs(i);
1983
1984 error_residual.norm = error_res.l2_norm();
1985 error_residual.u = error_res.block(u_dof).l2_norm();
1986 error_residual.p = error_res.block(p_dof).l2_norm();
1987 error_residual.J = error_res.block(J_dof).l2_norm();
1988 }
1989
1990
1991 // @sect4{Solid::get_error_update}
1992
1993 // Determine the true Newton update error for the problem
1994 template <int dim>
get_error_update(const BlockVector<double> & newton_update,Errors & error_update)1995 void Solid<dim>::get_error_update(const BlockVector<double> &newton_update,
1996 Errors & error_update)
1997 {
1998 BlockVector<double> error_ud(dofs_per_block);
1999 for (unsigned int i = 0; i < dof_handler.n_dofs(); ++i)
2000 if (!constraints.is_constrained(i))
2001 error_ud(i) = newton_update(i);
2002
2003 error_update.norm = error_ud.l2_norm();
2004 error_update.u = error_ud.block(u_dof).l2_norm();
2005 error_update.p = error_ud.block(p_dof).l2_norm();
2006 error_update.J = error_ud.block(J_dof).l2_norm();
2007 }
2008
2009
2010
2011 // @sect4{Solid::get_total_solution}
2012
2013 // This function provides the total solution, which is valid at any Newton
2014 // step. This is required as, to reduce computational error, the total
2015 // solution is only updated at the end of the timestep.
2016 template <int dim>
get_total_solution(const BlockVector<double> & solution_delta) const2017 BlockVector<double> Solid<dim>::get_total_solution(
2018 const BlockVector<double> &solution_delta) const
2019 {
2020 BlockVector<double> solution_total(solution_n);
2021 solution_total += solution_delta;
2022 return solution_total;
2023 }
2024
2025
2026 // @sect4{Solid::assemble_system_tangent}
2027
2028 // Since we use TBB for assembly, we simply setup a copy of the
2029 // data structures required for the process and pass them, along
2030 // with the memory addresses of the assembly functions to the
2031 // WorkStream object for processing. Note that we must ensure that
2032 // the matrix is reset before any assembly operations can occur.
2033 template <int dim>
assemble_system_tangent()2034 void Solid<dim>::assemble_system_tangent()
2035 {
2036 timer.enter_subsection("Assemble tangent matrix");
2037 std::cout << " ASM_K " << std::flush;
2038
2039 tangent_matrix = 0.0;
2040
2041 const UpdateFlags uf_cell(update_values | update_gradients |
2042 update_JxW_values);
2043
2044 PerTaskData_K per_task_data(dofs_per_cell);
2045 ScratchData_K scratch_data(fe, qf_cell, uf_cell);
2046
2047 // The syntax used here to pass data to the WorkStream class
2048 // is discussed in step-14. We need to use this particular
2049 // call to WorkStream because assemble_system_tangent_one_cell
2050 // is a constant function and copy_local_to_global_K is
2051 // non-constant.
2052 WorkStream::run(
2053 dof_handler.active_cell_iterators(),
2054 [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
2055 ScratchData_K & scratch,
2056 PerTaskData_K & data) {
2057 this->assemble_system_tangent_one_cell(cell, scratch, data);
2058 },
2059 [this](const PerTaskData_K &data) { this->copy_local_to_global_K(data); },
2060 scratch_data,
2061 per_task_data);
2062
2063 timer.leave_subsection();
2064 }
2065
2066 // This function adds the local contribution to the system matrix.
2067 // Note that we choose not to use the constraint matrix to do the
2068 // job for us because the tangent matrix and residual processes have
2069 // been split up into two separate functions.
2070 template <int dim>
copy_local_to_global_K(const PerTaskData_K & data)2071 void Solid<dim>::copy_local_to_global_K(const PerTaskData_K &data)
2072 {
2073 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2074 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2075 tangent_matrix.add(data.local_dof_indices[i],
2076 data.local_dof_indices[j],
2077 data.cell_matrix(i, j));
2078 }
2079
2080 // Of course, we still have to define how we assemble the tangent matrix
2081 // contribution for a single cell. We first need to reset and initialize some
2082 // of the scratch data structures and retrieve some basic information
2083 // regarding the DOF numbering on this cell. We can precalculate the cell
2084 // shape function values and gradients. Note that the shape function gradients
2085 // are defined with regard to the current configuration. That is
2086 // $\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi}
2087 // \ \mathbf{F}^{-1}$.
2088 template <int dim>
assemble_system_tangent_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,ScratchData_K & scratch,PerTaskData_K & data) const2089 void Solid<dim>::assemble_system_tangent_one_cell(
2090 const typename DoFHandler<dim>::active_cell_iterator &cell,
2091 ScratchData_K & scratch,
2092 PerTaskData_K & data) const
2093 {
2094 data.reset();
2095 scratch.reset();
2096 scratch.fe_values.reinit(cell);
2097 cell->get_dof_indices(data.local_dof_indices);
2098
2099 const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
2100 quadrature_point_history.get_data(cell);
2101 Assert(lqph.size() == n_q_points, ExcInternalError());
2102
2103 for (const unsigned int q_point :
2104 scratch.fe_values.quadrature_point_indices())
2105 {
2106 const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
2107 for (const unsigned int k : scratch.fe_values.dof_indices())
2108 {
2109 const unsigned int k_group = fe.system_to_base_index(k).first.first;
2110
2111 if (k_group == u_dof)
2112 {
2113 scratch.grad_Nx[q_point][k] =
2114 scratch.fe_values[u_fe].gradient(k, q_point) * F_inv;
2115 scratch.symm_grad_Nx[q_point][k] =
2116 symmetrize(scratch.grad_Nx[q_point][k]);
2117 }
2118 else if (k_group == p_dof)
2119 scratch.Nx[q_point][k] =
2120 scratch.fe_values[p_fe].value(k, q_point);
2121 else if (k_group == J_dof)
2122 scratch.Nx[q_point][k] =
2123 scratch.fe_values[J_fe].value(k, q_point);
2124 else
2125 Assert(k_group <= J_dof, ExcInternalError());
2126 }
2127 }
2128
2129 // Now we build the local cell stiffness matrix. Since the global and
2130 // local system matrices are symmetric, we can exploit this property by
2131 // building only the lower half of the local matrix and copying the values
2132 // to the upper half. So we only assemble half of the
2133 // $\mathsf{\mathbf{k}}_{uu}$, $\mathsf{\mathbf{k}}_{\widetilde{p}
2134 // \widetilde{p}} = \mathbf{0}$, $\mathsf{\mathbf{k}}_{\widetilde{J}
2135 // \widetilde{J}}$ blocks, while the whole
2136 // $\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}$,
2137 // $\mathsf{\mathbf{k}}_{u \widetilde{J}} = \mathbf{0}$,
2138 // $\mathsf{\mathbf{k}}_{u \widetilde{p}}$ blocks are built.
2139 //
2140 // In doing so, we first extract some configuration dependent variables
2141 // from our quadrature history objects for the current quadrature point.
2142 for (const unsigned int q_point :
2143 scratch.fe_values.quadrature_point_indices())
2144 {
2145 const Tensor<2, dim> tau = lqph[q_point]->get_tau();
2146 const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc();
2147 const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2();
2148 const double det_F = lqph[q_point]->get_det_F();
2149 const SymmetricTensor<2, dim> &I =
2150 Physics::Elasticity::StandardTensors<dim>::I;
2151
2152 // Next we define some aliases to make the assembly process easier to
2153 // follow
2154 const std::vector<double> & N = scratch.Nx[q_point];
2155 const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
2156 scratch.symm_grad_Nx[q_point];
2157 const std::vector<Tensor<2, dim>> &grad_Nx = scratch.grad_Nx[q_point];
2158 const double JxW = scratch.fe_values.JxW(q_point);
2159
2160 for (const unsigned int i : scratch.fe_values.dof_indices())
2161 {
2162 const unsigned int component_i =
2163 fe.system_to_component_index(i).first;
2164 const unsigned int i_group = fe.system_to_base_index(i).first.first;
2165
2166 for (const unsigned int j :
2167 scratch.fe_values.dof_indices_ending_at(i))
2168 {
2169 const unsigned int component_j =
2170 fe.system_to_component_index(j).first;
2171 const unsigned int j_group =
2172 fe.system_to_base_index(j).first.first;
2173
2174 // This is the $\mathsf{\mathbf{k}}_{uu}$
2175 // contribution. It comprises a material contribution, and a
2176 // geometrical stress contribution which is only added along
2177 // the local matrix diagonals:
2178 if ((i_group == j_group) && (i_group == u_dof))
2179 {
2180 // The material contribution:
2181 data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc * //
2182 symm_grad_Nx[j] * JxW; //
2183
2184 // The geometrical stress contribution:
2185 if (component_i == component_j)
2186 data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau *
2187 grad_Nx[j][component_j] * JxW;
2188 }
2189 // Next is the $\mathsf{\mathbf{k}}_{ \widetilde{p} u}$
2190 // contribution
2191 else if ((i_group == p_dof) && (j_group == u_dof))
2192 {
2193 data.cell_matrix(i, j) += N[i] * det_F * //
2194 (symm_grad_Nx[j] * I) * JxW; //
2195 }
2196 // and lastly the $\mathsf{\mathbf{k}}_{ \widetilde{J}
2197 // \widetilde{p}}$ and $\mathsf{\mathbf{k}}_{ \widetilde{J}
2198 // \widetilde{J}}$ contributions:
2199 else if ((i_group == J_dof) && (j_group == p_dof))
2200 data.cell_matrix(i, j) -= N[i] * N[j] * JxW;
2201 else if ((i_group == j_group) && (i_group == J_dof))
2202 data.cell_matrix(i, j) += N[i] * d2Psi_vol_dJ2 * N[j] * JxW;
2203 else
2204 Assert((i_group <= J_dof) && (j_group <= J_dof),
2205 ExcInternalError());
2206 }
2207 }
2208 }
2209
2210 // Finally, we need to copy the lower half of the local matrix into the
2211 // upper half:
2212 for (const unsigned int i : scratch.fe_values.dof_indices())
2213 for (const unsigned int j :
2214 scratch.fe_values.dof_indices_starting_at(i + 1))
2215 data.cell_matrix(i, j) = data.cell_matrix(j, i);
2216 }
2217
2218 // @sect4{Solid::assemble_system_rhs}
2219 // The assembly of the right-hand side process is similar to the
2220 // tangent matrix, so we will not describe it in too much detail.
2221 // Note that since we are describing a problem with Neumann BCs,
2222 // we will need the face normals and so must specify this in the
2223 // update flags.
2224 template <int dim>
assemble_system_rhs()2225 void Solid<dim>::assemble_system_rhs()
2226 {
2227 timer.enter_subsection("Assemble system right-hand side");
2228 std::cout << " ASM_R " << std::flush;
2229
2230 system_rhs = 0.0;
2231
2232 const UpdateFlags uf_cell(update_values | update_gradients |
2233 update_JxW_values);
2234 const UpdateFlags uf_face(update_values | update_normal_vectors |
2235 update_JxW_values);
2236
2237 PerTaskData_RHS per_task_data(dofs_per_cell);
2238 ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face);
2239
2240 WorkStream::run(
2241 dof_handler.active_cell_iterators(),
2242 [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
2243 ScratchData_RHS & scratch,
2244 PerTaskData_RHS & data) {
2245 this->assemble_system_rhs_one_cell(cell, scratch, data);
2246 },
2247 [this](const PerTaskData_RHS &data) {
2248 this->copy_local_to_global_rhs(data);
2249 },
2250 scratch_data,
2251 per_task_data);
2252
2253 timer.leave_subsection();
2254 }
2255
2256
2257
2258 template <int dim>
copy_local_to_global_rhs(const PerTaskData_RHS & data)2259 void Solid<dim>::copy_local_to_global_rhs(const PerTaskData_RHS &data)
2260 {
2261 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2262 system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i);
2263 }
2264
2265
2266
2267 template <int dim>
assemble_system_rhs_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,ScratchData_RHS & scratch,PerTaskData_RHS & data) const2268 void Solid<dim>::assemble_system_rhs_one_cell(
2269 const typename DoFHandler<dim>::active_cell_iterator &cell,
2270 ScratchData_RHS & scratch,
2271 PerTaskData_RHS & data) const
2272 {
2273 data.reset();
2274 scratch.reset();
2275 scratch.fe_values.reinit(cell);
2276 cell->get_dof_indices(data.local_dof_indices);
2277
2278 const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
2279 quadrature_point_history.get_data(cell);
2280 Assert(lqph.size() == n_q_points, ExcInternalError());
2281
2282 for (const unsigned int q_point :
2283 scratch.fe_values.quadrature_point_indices())
2284 {
2285 const Tensor<2, dim> F_inv = lqph[q_point]->get_F_inv();
2286
2287 for (const unsigned int k : scratch.fe_values.dof_indices())
2288 {
2289 const unsigned int k_group = fe.system_to_base_index(k).first.first;
2290
2291 if (k_group == u_dof)
2292 scratch.symm_grad_Nx[q_point][k] = symmetrize(
2293 scratch.fe_values[u_fe].gradient(k, q_point) * F_inv);
2294 else if (k_group == p_dof)
2295 scratch.Nx[q_point][k] =
2296 scratch.fe_values[p_fe].value(k, q_point);
2297 else if (k_group == J_dof)
2298 scratch.Nx[q_point][k] =
2299 scratch.fe_values[J_fe].value(k, q_point);
2300 else
2301 Assert(k_group <= J_dof, ExcInternalError());
2302 }
2303 }
2304
2305 for (const unsigned int q_point :
2306 scratch.fe_values.quadrature_point_indices())
2307 {
2308 const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau();
2309 const double det_F = lqph[q_point]->get_det_F();
2310 const double J_tilde = lqph[q_point]->get_J_tilde();
2311 const double p_tilde = lqph[q_point]->get_p_tilde();
2312 const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ();
2313
2314 const std::vector<double> & N = scratch.Nx[q_point];
2315 const std::vector<SymmetricTensor<2, dim>> &symm_grad_Nx =
2316 scratch.symm_grad_Nx[q_point];
2317 const double JxW = scratch.fe_values.JxW(q_point);
2318
2319 // We first compute the contributions
2320 // from the internal forces. Note, by
2321 // definition of the rhs as the negative
2322 // of the residual, these contributions
2323 // are subtracted.
2324 for (const unsigned int i : scratch.fe_values.dof_indices())
2325 {
2326 const unsigned int i_group = fe.system_to_base_index(i).first.first;
2327
2328 if (i_group == u_dof)
2329 data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW;
2330 else if (i_group == p_dof)
2331 data.cell_rhs(i) -= N[i] * (det_F - J_tilde) * JxW;
2332 else if (i_group == J_dof)
2333 data.cell_rhs(i) -= N[i] * (dPsi_vol_dJ - p_tilde) * JxW;
2334 else
2335 Assert(i_group <= J_dof, ExcInternalError());
2336 }
2337 }
2338
2339 // Next we assemble the Neumann contribution. We first check to see it the
2340 // cell face exists on a boundary on which a traction is applied and add
2341 // the contribution if this is the case.
2342 for (const auto &face : cell->face_iterators())
2343 if (face->at_boundary() && face->boundary_id() == 6)
2344 {
2345 scratch.fe_face_values.reinit(cell, face);
2346
2347 for (const unsigned int f_q_point :
2348 scratch.fe_face_values.quadrature_point_indices())
2349 {
2350 const Tensor<1, dim> &N =
2351 scratch.fe_face_values.normal_vector(f_q_point);
2352
2353 // Using the face normal at this quadrature point we specify the
2354 // traction in reference configuration. For this problem, a
2355 // defined pressure is applied in the reference configuration.
2356 // The direction of the applied traction is assumed not to
2357 // evolve with the deformation of the domain. The traction is
2358 // defined using the first Piola-Kirchhoff stress is simply
2359 // $\mathbf{t} = \mathbf{P}\mathbf{N} = [p_0 \mathbf{I}]
2360 // \mathbf{N} = p_0 \mathbf{N}$ We use the time variable to
2361 // linearly ramp up the pressure load.
2362 //
2363 // Note that the contributions to the right hand side vector we
2364 // compute here only exist in the displacement components of the
2365 // vector.
2366 static const double p0 =
2367 -4.0 / (parameters.scale * parameters.scale);
2368 const double time_ramp = (time.current() / time.end());
2369 const double pressure = p0 * parameters.p_p0 * time_ramp;
2370 const Tensor<1, dim> traction = pressure * N;
2371
2372 for (const unsigned int i : scratch.fe_values.dof_indices())
2373 {
2374 const unsigned int i_group =
2375 fe.system_to_base_index(i).first.first;
2376
2377 if (i_group == u_dof)
2378 {
2379 const unsigned int component_i =
2380 fe.system_to_component_index(i).first;
2381 const double Ni =
2382 scratch.fe_face_values.shape_value(i, f_q_point);
2383 const double JxW = scratch.fe_face_values.JxW(f_q_point);
2384
2385 data.cell_rhs(i) += (Ni * traction[component_i]) * JxW;
2386 }
2387 }
2388 }
2389 }
2390 }
2391
2392 // @sect4{Solid::make_constraints}
2393 // The constraints for this problem are simple to describe.
2394 // However, since we are dealing with an iterative Newton method,
2395 // it should be noted that any displacement constraints should only
2396 // be specified at the zeroth iteration and subsequently no
2397 // additional contributions are to be made since the constraints
2398 // are already exactly satisfied.
2399 template <int dim>
make_constraints(const int & it_nr)2400 void Solid<dim>::make_constraints(const int &it_nr)
2401 {
2402 std::cout << " CST " << std::flush;
2403
2404 // Since the constraints are different at different Newton iterations, we
2405 // need to clear the constraints matrix and completely rebuild
2406 // it. However, after the first iteration, the constraints remain the same
2407 // and we can simply skip the rebuilding step if we do not clear it.
2408 if (it_nr > 1)
2409 return;
2410 constraints.clear();
2411 const bool apply_dirichlet_bc = (it_nr == 0);
2412
2413 // In this particular example, the boundary values will be calculated for
2414 // the two first iterations of Newton's algorithm. In general, one would
2415 // build non-homogeneous constraints in the zeroth iteration (that is, when
2416 // `apply_dirichlet_bc == true`) and build only the corresponding
2417 // homogeneous constraints in the following step. While the current
2418 // example has only homogeneous constraints, previous experiences have
2419 // shown that a common error is forgetting to add the extra condition when
2420 // refactoring the code to specific uses. This could lead errors that are
2421 // hard to debug. In this spirit, we choose to make the code more verbose
2422 // in terms of what operations are performed at each Newton step.
2423 //
2424 // The boundary conditions for the indentation problem are as follows: On
2425 // the -x, -y and -z faces (IDs 0,2,4) we set up a symmetry condition to
2426 // allow only planar movement while the +x and +z faces (IDs 1,5) are
2427 // traction free. In this contrived problem, part of the +y face (ID 3) is
2428 // set to have no motion in the x- and z-component. Finally, as described
2429 // earlier, the other part of the +y face has an the applied pressure but
2430 // is also constrained in the x- and z-directions.
2431 //
2432 // In the following, we will have to tell the function interpolation
2433 // boundary values which components of the solution vector should be
2434 // constrained (i.e., whether it's the x-, y-, z-displacements or
2435 // combinations thereof). This is done using ComponentMask objects (see
2436 // @ref GlossComponentMask) which we can get from the finite element if we
2437 // provide it with an extractor object for the component we wish to
2438 // select. To this end we first set up such extractor objects and later
2439 // use it when generating the relevant component masks:
2440 const FEValuesExtractors::Scalar x_displacement(0);
2441 const FEValuesExtractors::Scalar y_displacement(1);
2442
2443 {
2444 const int boundary_id = 0;
2445
2446 if (apply_dirichlet_bc == true)
2447 VectorTools::interpolate_boundary_values(
2448 dof_handler,
2449 boundary_id,
2450 Functions::ZeroFunction<dim>(n_components),
2451 constraints,
2452 fe.component_mask(x_displacement));
2453 else
2454 VectorTools::interpolate_boundary_values(
2455 dof_handler,
2456 boundary_id,
2457 Functions::ZeroFunction<dim>(n_components),
2458 constraints,
2459 fe.component_mask(x_displacement));
2460 }
2461 {
2462 const int boundary_id = 2;
2463
2464 if (apply_dirichlet_bc == true)
2465 VectorTools::interpolate_boundary_values(
2466 dof_handler,
2467 boundary_id,
2468 Functions::ZeroFunction<dim>(n_components),
2469 constraints,
2470 fe.component_mask(y_displacement));
2471 else
2472 VectorTools::interpolate_boundary_values(
2473 dof_handler,
2474 boundary_id,
2475 Functions::ZeroFunction<dim>(n_components),
2476 constraints,
2477 fe.component_mask(y_displacement));
2478 }
2479
2480 if (dim == 3)
2481 {
2482 const FEValuesExtractors::Scalar z_displacement(2);
2483
2484 {
2485 const int boundary_id = 3;
2486
2487 if (apply_dirichlet_bc == true)
2488 VectorTools::interpolate_boundary_values(
2489 dof_handler,
2490 boundary_id,
2491 Functions::ZeroFunction<dim>(n_components),
2492 constraints,
2493 (fe.component_mask(x_displacement) |
2494 fe.component_mask(z_displacement)));
2495 else
2496 VectorTools::interpolate_boundary_values(
2497 dof_handler,
2498 boundary_id,
2499 Functions::ZeroFunction<dim>(n_components),
2500 constraints,
2501 (fe.component_mask(x_displacement) |
2502 fe.component_mask(z_displacement)));
2503 }
2504 {
2505 const int boundary_id = 4;
2506
2507 if (apply_dirichlet_bc == true)
2508 VectorTools::interpolate_boundary_values(
2509 dof_handler,
2510 boundary_id,
2511 Functions::ZeroFunction<dim>(n_components),
2512 constraints,
2513 fe.component_mask(z_displacement));
2514 else
2515 VectorTools::interpolate_boundary_values(
2516 dof_handler,
2517 boundary_id,
2518 Functions::ZeroFunction<dim>(n_components),
2519 constraints,
2520 fe.component_mask(z_displacement));
2521 }
2522
2523 {
2524 const int boundary_id = 6;
2525
2526 if (apply_dirichlet_bc == true)
2527 VectorTools::interpolate_boundary_values(
2528 dof_handler,
2529 boundary_id,
2530 Functions::ZeroFunction<dim>(n_components),
2531 constraints,
2532 (fe.component_mask(x_displacement) |
2533 fe.component_mask(z_displacement)));
2534 else
2535 VectorTools::interpolate_boundary_values(
2536 dof_handler,
2537 boundary_id,
2538 Functions::ZeroFunction<dim>(n_components),
2539 constraints,
2540 (fe.component_mask(x_displacement) |
2541 fe.component_mask(z_displacement)));
2542 }
2543 }
2544 else
2545 {
2546 {
2547 const int boundary_id = 3;
2548
2549 if (apply_dirichlet_bc == true)
2550 VectorTools::interpolate_boundary_values(
2551 dof_handler,
2552 boundary_id,
2553 Functions::ZeroFunction<dim>(n_components),
2554 constraints,
2555 (fe.component_mask(x_displacement)));
2556 else
2557 VectorTools::interpolate_boundary_values(
2558 dof_handler,
2559 boundary_id,
2560 Functions::ZeroFunction<dim>(n_components),
2561 constraints,
2562 (fe.component_mask(x_displacement)));
2563 }
2564 {
2565 const int boundary_id = 6;
2566
2567 if (apply_dirichlet_bc == true)
2568 VectorTools::interpolate_boundary_values(
2569 dof_handler,
2570 boundary_id,
2571 Functions::ZeroFunction<dim>(n_components),
2572 constraints,
2573 (fe.component_mask(x_displacement)));
2574 else
2575 VectorTools::interpolate_boundary_values(
2576 dof_handler,
2577 boundary_id,
2578 Functions::ZeroFunction<dim>(n_components),
2579 constraints,
2580 (fe.component_mask(x_displacement)));
2581 }
2582 }
2583
2584 constraints.close();
2585 }
2586
2587 // @sect4{Solid::assemble_sc}
2588 // Solving the entire block system is a bit problematic as there are no
2589 // contributions to the $\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}$
2590 // block, rendering it noninvertible (when using an iterative solver).
2591 // Since the pressure and dilatation variables DOFs are discontinuous, we can
2592 // condense them out to form a smaller displacement-only system which
2593 // we will then solve and subsequently post-process to retrieve the
2594 // pressure and dilatation solutions.
2595
2596 // The static condensation process could be performed at a global level but we
2597 // need the inverse of one of the blocks. However, since the pressure and
2598 // dilatation variables are discontinuous, the static condensation (SC)
2599 // operation can also be done on a per-cell basis and we can produce the
2600 // inverse of the block-diagonal
2601 // $\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}$ block by inverting the
2602 // local blocks. We can again use TBB to do this since each operation will be
2603 // independent of one another.
2604 //
2605 // Using the TBB via the WorkStream class, we assemble the contributions to
2606 // form
2607 // $
2608 // \mathsf{\mathbf{K}}_{\textrm{con}}
2609 // = \bigl[ \mathsf{\mathbf{K}}_{uu} +
2610 // \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
2611 // $
2612 // from each element's contributions. These contributions are then added to
2613 // the global stiffness matrix. Given this description, the following two
2614 // functions should be clear:
2615 template <int dim>
assemble_sc()2616 void Solid<dim>::assemble_sc()
2617 {
2618 timer.enter_subsection("Perform static condensation");
2619 std::cout << " ASM_SC " << std::flush;
2620
2621 PerTaskData_SC per_task_data(dofs_per_cell,
2622 element_indices_u.size(),
2623 element_indices_p.size(),
2624 element_indices_J.size());
2625 ScratchData_SC scratch_data;
2626
2627 WorkStream::run(dof_handler.active_cell_iterators(),
2628 *this,
2629 &Solid::assemble_sc_one_cell,
2630 &Solid::copy_local_to_global_sc,
2631 scratch_data,
2632 per_task_data);
2633
2634 timer.leave_subsection();
2635 }
2636
2637
2638 template <int dim>
copy_local_to_global_sc(const PerTaskData_SC & data)2639 void Solid<dim>::copy_local_to_global_sc(const PerTaskData_SC &data)
2640 {
2641 for (unsigned int i = 0; i < dofs_per_cell; ++i)
2642 for (unsigned int j = 0; j < dofs_per_cell; ++j)
2643 tangent_matrix.add(data.local_dof_indices[i],
2644 data.local_dof_indices[j],
2645 data.cell_matrix(i, j));
2646 }
2647
2648
2649 // Now we describe the static condensation process. As per usual, we must
2650 // first find out which global numbers the degrees of freedom on this cell
2651 // have and reset some data structures:
2652 template <int dim>
assemble_sc_one_cell(const typename DoFHandler<dim>::active_cell_iterator & cell,ScratchData_SC & scratch,PerTaskData_SC & data)2653 void Solid<dim>::assemble_sc_one_cell(
2654 const typename DoFHandler<dim>::active_cell_iterator &cell,
2655 ScratchData_SC & scratch,
2656 PerTaskData_SC & data)
2657 {
2658 data.reset();
2659 scratch.reset();
2660 cell->get_dof_indices(data.local_dof_indices);
2661
2662 // We now extract the contribution of the dofs associated with the current
2663 // cell to the global stiffness matrix. The discontinuous nature of the
2664 // $\widetilde{p}$ and $\widetilde{J}$ interpolations mean that their is
2665 // no coupling of the local contributions at the global level. This is not
2666 // the case with the $\mathbf{u}$ dof. In other words,
2667 // $\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}$,
2668 // $\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{p}}$ and
2669 // $\mathsf{\mathbf{k}}_{\widetilde{J} \widetilde{p}}$, when extracted
2670 // from the global stiffness matrix are the element contributions. This
2671 // is not the case for $\mathsf{\mathbf{k}}_{uu}$.
2672 //
2673 // Note: A lower-case symbol is used to denote element stiffness matrices.
2674
2675 // Currently the matrix corresponding to
2676 // the dof associated with the current element
2677 // (denoted somewhat loosely as $\mathsf{\mathbf{k}}$)
2678 // is of the form:
2679 // @f{align*}
2680 // \begin{bmatrix}
2681 // \mathsf{\mathbf{k}}_{uu} & \mathsf{\mathbf{k}}_{u\widetilde{p}}
2682 // & \mathbf{0}
2683 // \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
2684 // \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}
2685 // \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
2686 // \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
2687 // @f}
2688 //
2689 // We now need to modify it such that it appear as
2690 // @f{align*}
2691 // \begin{bmatrix}
2692 // \mathsf{\mathbf{k}}_{\textrm{con}} &
2693 // \mathsf{\mathbf{k}}_{u\widetilde{p}} & \mathbf{0}
2694 // \\ \mathsf{\mathbf{k}}_{\widetilde{p}u} & \mathbf{0} &
2695 // \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
2696 // \\ \mathbf{0} & \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}} &
2697 // \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}} \end{bmatrix}
2698 // @f}
2699 // with $\mathsf{\mathbf{k}}_{\textrm{con}} = \bigl[
2700 // \mathsf{\mathbf{k}}_{uu} +\overline{\overline{\mathsf{\mathbf{k}}}}~
2701 // \bigr]$ where $ \overline{\overline{\mathsf{\mathbf{k}}}}
2702 // \dealcoloneq \mathsf{\mathbf{k}}_{u\widetilde{p}}
2703 // \overline{\mathsf{\mathbf{k}}} \mathsf{\mathbf{k}}_{\widetilde{p}u}
2704 // $
2705 // and
2706 // $
2707 // \overline{\mathsf{\mathbf{k}}} =
2708 // \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{p}}^{-1}
2709 // \mathsf{\mathbf{k}}_{\widetilde{J}\widetilde{J}}
2710 // \mathsf{\mathbf{k}}_{\widetilde{p}\widetilde{J}}^{-1}
2711 // $.
2712 //
2713 // At this point, we need to take note of
2714 // the fact that global data already exists
2715 // in the $\mathsf{\mathbf{K}}_{uu}$,
2716 // $\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}$
2717 // and
2718 // $\mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{p}}$
2719 // sub-blocks. So if we are to modify them, we must account for the data
2720 // that is already there (i.e. simply add to it or remove it if
2721 // necessary). Since the copy_local_to_global operation is a "+="
2722 // operation, we need to take this into account
2723 //
2724 // For the $\mathsf{\mathbf{K}}_{uu}$ block in particular, this means that
2725 // contributions have been added from the surrounding cells, so we need to
2726 // be careful when we manipulate this block. We can't just erase the
2727 // sub-blocks.
2728 //
2729 // This is the strategy we will employ to get the sub-blocks we want:
2730 //
2731 // - $ {\mathsf{\mathbf{k}}}_{\textrm{store}}$:
2732 // Since we don't have access to $\mathsf{\mathbf{k}}_{uu}$,
2733 // but we know its contribution is added to
2734 // the global $\mathsf{\mathbf{K}}_{uu}$ matrix, we just want
2735 // to add the element wise
2736 // static-condensation $\overline{\overline{\mathsf{\mathbf{k}}}}$.
2737 //
2738 // - $\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}$:
2739 // Similarly, $\mathsf{\mathbf{k}}_{\widetilde{p}
2740 // \widetilde{J}}$ exists in
2741 // the subblock. Since the copy
2742 // operation is a += operation, we
2743 // need to subtract the existing
2744 // $\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}$
2745 // submatrix in addition to
2746 // "adding" that which we wish to
2747 // replace it with.
2748 //
2749 // - $\mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}$:
2750 // Since the global matrix
2751 // is symmetric, this block is the
2752 // same as the one above and we
2753 // can simply use
2754 // $\mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}$
2755 // as a substitute for this one.
2756 //
2757 // We first extract element data from the
2758 // system matrix. So first we get the
2759 // entire subblock for the cell, then
2760 // extract $\mathsf{\mathbf{k}}$
2761 // for the dofs associated with
2762 // the current element
2763 data.k_orig.extract_submatrix_from(tangent_matrix,
2764 data.local_dof_indices,
2765 data.local_dof_indices);
2766 // and next the local matrices for
2767 // $\mathsf{\mathbf{k}}_{ \widetilde{p} u}$
2768 // $\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}$
2769 // and
2770 // $\mathsf{\mathbf{k}}_{ \widetilde{J} \widetilde{J}}$:
2771 data.k_pu.extract_submatrix_from(data.k_orig,
2772 element_indices_p,
2773 element_indices_u);
2774 data.k_pJ.extract_submatrix_from(data.k_orig,
2775 element_indices_p,
2776 element_indices_J);
2777 data.k_JJ.extract_submatrix_from(data.k_orig,
2778 element_indices_J,
2779 element_indices_J);
2780
2781 // To get the inverse of $\mathsf{\mathbf{k}}_{\widetilde{p}
2782 // \widetilde{J}}$, we invert it directly. This operation is relatively
2783 // inexpensive since $\mathsf{\mathbf{k}}_{\widetilde{p} \widetilde{J}}$
2784 // since block-diagonal.
2785 data.k_pJ_inv.invert(data.k_pJ);
2786
2787 // Now we can make condensation terms to
2788 // add to the $\mathsf{\mathbf{k}}_{uu}$
2789 // block and put them in
2790 // the cell local matrix
2791 // $
2792 // \mathsf{\mathbf{A}}
2793 // =
2794 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
2795 // \mathsf{\mathbf{k}}_{\widetilde{p} u}
2796 // $:
2797 data.k_pJ_inv.mmult(data.A, data.k_pu);
2798 // $
2799 // \mathsf{\mathbf{B}}
2800 // =
2801 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
2802 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
2803 // \mathsf{\mathbf{k}}_{\widetilde{p} u}
2804 // $
2805 data.k_JJ.mmult(data.B, data.A);
2806 // $
2807 // \mathsf{\mathbf{C}}
2808 // =
2809 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
2810 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
2811 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
2812 // \mathsf{\mathbf{k}}_{\widetilde{p} u}
2813 // $
2814 data.k_pJ_inv.Tmmult(data.C, data.B);
2815 // $
2816 // \overline{\overline{\mathsf{\mathbf{k}}}}
2817 // =
2818 // \mathsf{\mathbf{k}}_{u \widetilde{p}}
2819 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{p}}
2820 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{J} \widetilde{J}}
2821 // \mathsf{\mathbf{k}}^{-1}_{\widetilde{p} \widetilde{J}}
2822 // \mathsf{\mathbf{k}}_{\widetilde{p} u}
2823 // $
2824 data.k_pu.Tmmult(data.k_bbar, data.C);
2825 data.k_bbar.scatter_matrix_to(element_indices_u,
2826 element_indices_u,
2827 data.cell_matrix);
2828
2829 // Next we place
2830 // $\mathsf{\mathbf{k}}^{-1}_{ \widetilde{p} \widetilde{J}}$
2831 // in the
2832 // $\mathsf{\mathbf{k}}_{ \widetilde{p} \widetilde{J}}$
2833 // block for post-processing. Note again
2834 // that we need to remove the
2835 // contribution that already exists there.
2836 data.k_pJ_inv.add(-1.0, data.k_pJ);
2837 data.k_pJ_inv.scatter_matrix_to(element_indices_p,
2838 element_indices_J,
2839 data.cell_matrix);
2840 }
2841
2842 // @sect4{Solid::solve_linear_system}
2843 // We now have all of the necessary components to use one of two possible
2844 // methods to solve the linearised system. The first is to perform static
2845 // condensation on an element level, which requires some alterations
2846 // to the tangent matrix and RHS vector. Alternatively, the full block
2847 // system can be solved by performing condensation on a global level.
2848 // Below we implement both approaches.
2849 template <int dim>
2850 std::pair<unsigned int, double>
solve_linear_system(BlockVector<double> & newton_update)2851 Solid<dim>::solve_linear_system(BlockVector<double> &newton_update)
2852 {
2853 unsigned int lin_it = 0;
2854 double lin_res = 0.0;
2855
2856 if (parameters.use_static_condensation == true)
2857 {
2858 // Firstly, here is the approach using the (permanent) augmentation of
2859 // the tangent matrix. For the following, recall that
2860 // @f{align*}
2861 // \mathsf{\mathbf{K}}_{\textrm{store}}
2862 //\dealcoloneq
2863 // \begin{bmatrix}
2864 // \mathsf{\mathbf{K}}_{\textrm{con}} &
2865 // \mathsf{\mathbf{K}}_{u\widetilde{p}} & \mathbf{0}
2866 // \\ \mathsf{\mathbf{K}}_{\widetilde{p}u} & \mathbf{0} &
2867 // \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
2868 // \\ \mathbf{0} &
2869 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}} &
2870 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}} \end{bmatrix} \, .
2871 // @f}
2872 // and
2873 // @f{align*}
2874 // d \widetilde{\mathsf{\mathbf{p}}}
2875 // & =
2876 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
2877 // \bigl[
2878 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2879 // -
2880 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
2881 // d \widetilde{\mathsf{\mathbf{J}}} \bigr]
2882 // \\ d \widetilde{\mathsf{\mathbf{J}}}
2883 // & =
2884 // \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}
2885 // \bigl[
2886 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2887 // - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
2888 // \mathsf{\mathbf{u}} \bigr]
2889 // \\ \Rightarrow d \widetilde{\mathsf{\mathbf{p}}}
2890 // &= \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
2891 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2892 // -
2893 // \underbrace{\bigl[\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
2894 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
2895 // \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}\bigr]}_{\overline{\mathsf{\mathbf{K}}}}\bigl[
2896 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2897 // - \mathsf{\mathbf{K}}_{\widetilde{p}u} d
2898 // \mathsf{\mathbf{u}} \bigr]
2899 // @f}
2900 // and thus
2901 // @f[
2902 // \underbrace{\bigl[ \mathsf{\mathbf{K}}_{uu} +
2903 // \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]
2904 // }_{\mathsf{\mathbf{K}}_{\textrm{con}}} d
2905 // \mathsf{\mathbf{u}}
2906 // =
2907 // \underbrace{
2908 // \Bigl[
2909 // \mathsf{\mathbf{F}}_{u}
2910 // - \mathsf{\mathbf{K}}_{u\widetilde{p}} \bigl[
2911 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
2912 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2913 // -
2914 // \overline{\mathsf{\mathbf{K}}}\mathsf{\mathbf{F}}_{\widetilde{p}}
2915 // \bigr]
2916 // \Bigr]}_{\mathsf{\mathbf{F}}_{\textrm{con}}}
2917 // @f]
2918 // where
2919 // @f[
2920 // \overline{\overline{\mathsf{\mathbf{K}}}} \dealcoloneq
2921 // \mathsf{\mathbf{K}}_{u\widetilde{p}}
2922 // \overline{\mathsf{\mathbf{K}}}
2923 // \mathsf{\mathbf{K}}_{\widetilde{p}u} \, .
2924 // @f]
2925
2926 // At the top, we allocate two temporary vectors to help with the
2927 // static condensation, and variables to store the number of
2928 // linear solver iterations and the (hopefully converged) residual.
2929 BlockVector<double> A(dofs_per_block);
2930 BlockVector<double> B(dofs_per_block);
2931
2932
2933 // In the first step of this function, we solve for the incremental
2934 // displacement $d\mathbf{u}$. To this end, we perform static
2935 // condensation to make
2936 // $\mathsf{\mathbf{K}}_{\textrm{con}}
2937 // = \bigl[ \mathsf{\mathbf{K}}_{uu} +
2938 // \overline{\overline{\mathsf{\mathbf{K}}}}~ \bigr]$
2939 // and put
2940 // $\mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}$
2941 // in the original $\mathsf{\mathbf{K}}_{\widetilde{p} \widetilde{J}}$
2942 // block. That is, we make $\mathsf{\mathbf{K}}_{\textrm{store}}$.
2943 {
2944 assemble_sc();
2945
2946 // $
2947 // \mathsf{\mathbf{A}}_{\widetilde{J}}
2948 // =
2949 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
2950 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2951 // $
2952 tangent_matrix.block(p_dof, J_dof)
2953 .vmult(A.block(J_dof), system_rhs.block(p_dof));
2954 // $
2955 // \mathsf{\mathbf{B}}_{\widetilde{J}}
2956 // =
2957 // \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
2958 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
2959 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2960 // $
2961 tangent_matrix.block(J_dof, J_dof)
2962 .vmult(B.block(J_dof), A.block(J_dof));
2963 // $
2964 // \mathsf{\mathbf{A}}_{\widetilde{J}}
2965 // =
2966 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2967 // -
2968 // \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
2969 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
2970 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2971 // $
2972 A.block(J_dof) = system_rhs.block(J_dof);
2973 A.block(J_dof) -= B.block(J_dof);
2974 // $
2975 // \mathsf{\mathbf{A}}_{\widetilde{J}}
2976 // =
2977 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
2978 // [
2979 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2980 // -
2981 // \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
2982 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
2983 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2984 // ]
2985 // $
2986 tangent_matrix.block(p_dof, J_dof)
2987 .Tvmult(A.block(p_dof), A.block(J_dof));
2988 // $
2989 // \mathsf{\mathbf{A}}_{u}
2990 // =
2991 // \mathsf{\mathbf{K}}_{u \widetilde{p}}
2992 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
2993 // [
2994 // \mathsf{\mathbf{F}}_{\widetilde{J}}
2995 // -
2996 // \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
2997 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
2998 // \mathsf{\mathbf{F}}_{\widetilde{p}}
2999 // ]
3000 // $
3001 tangent_matrix.block(u_dof, p_dof)
3002 .vmult(A.block(u_dof), A.block(p_dof));
3003 // $
3004 // \mathsf{\mathbf{F}}_{\text{con}}
3005 // =
3006 // \mathsf{\mathbf{F}}_{u}
3007 // -
3008 // \mathsf{\mathbf{K}}_{u \widetilde{p}}
3009 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{J} \widetilde{p}}
3010 // [
3011 // \mathsf{\mathbf{F}}_{\widetilde{J}}
3012 // -
3013 // \mathsf{\mathbf{K}}_{\widetilde{J} \widetilde{J}}
3014 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p} \widetilde{J}}
3015 // \mathsf{\mathbf{F}}_{\widetilde{p}}
3016 // ]
3017 // $
3018 system_rhs.block(u_dof) -= A.block(u_dof);
3019
3020 timer.enter_subsection("Linear solver");
3021 std::cout << " SLV " << std::flush;
3022 if (parameters.type_lin == "CG")
3023 {
3024 const auto solver_its = static_cast<unsigned int>(
3025 tangent_matrix.block(u_dof, u_dof).m() *
3026 parameters.max_iterations_lin);
3027 const double tol_sol =
3028 parameters.tol_lin * system_rhs.block(u_dof).l2_norm();
3029
3030 SolverControl solver_control(solver_its, tol_sol);
3031
3032 GrowingVectorMemory<Vector<double>> GVM;
3033 SolverCG<Vector<double>> solver_CG(solver_control, GVM);
3034
3035 // We've chosen by default a SSOR preconditioner as it appears to
3036 // provide the fastest solver convergence characteristics for this
3037 // problem on a single-thread machine. However, this might not be
3038 // true for different problem sizes.
3039 PreconditionSelector<SparseMatrix<double>, Vector<double>>
3040 preconditioner(parameters.preconditioner_type,
3041 parameters.preconditioner_relaxation);
3042 preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof));
3043
3044 solver_CG.solve(tangent_matrix.block(u_dof, u_dof),
3045 newton_update.block(u_dof),
3046 system_rhs.block(u_dof),
3047 preconditioner);
3048
3049 lin_it = solver_control.last_step();
3050 lin_res = solver_control.last_value();
3051 }
3052 else if (parameters.type_lin == "Direct")
3053 {
3054 // Otherwise if the problem is small
3055 // enough, a direct solver can be
3056 // utilised.
3057 SparseDirectUMFPACK A_direct;
3058 A_direct.initialize(tangent_matrix.block(u_dof, u_dof));
3059 A_direct.vmult(newton_update.block(u_dof),
3060 system_rhs.block(u_dof));
3061
3062 lin_it = 1;
3063 lin_res = 0.0;
3064 }
3065 else
3066 Assert(false, ExcMessage("Linear solver type not implemented"));
3067
3068 timer.leave_subsection();
3069 }
3070
3071 // Now that we have the displacement update, distribute the constraints
3072 // back to the Newton update:
3073 constraints.distribute(newton_update);
3074
3075 timer.enter_subsection("Linear solver postprocessing");
3076 std::cout << " PP " << std::flush;
3077
3078 // The next step after solving the displacement
3079 // problem is to post-process to get the
3080 // dilatation solution from the
3081 // substitution:
3082 // $
3083 // d \widetilde{\mathsf{\mathbf{J}}}
3084 // = \mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1} \bigl[
3085 // \mathsf{\mathbf{F}}_{\widetilde{p}}
3086 // - \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
3087 // \bigr]
3088 // $
3089 {
3090 // $
3091 // \mathsf{\mathbf{A}}_{\widetilde{p}}
3092 // =
3093 // \mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
3094 // $
3095 tangent_matrix.block(p_dof, u_dof)
3096 .vmult(A.block(p_dof), newton_update.block(u_dof));
3097 // $
3098 // \mathsf{\mathbf{A}}_{\widetilde{p}}
3099 // =
3100 // -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
3101 // $
3102 A.block(p_dof) *= -1.0;
3103 // $
3104 // \mathsf{\mathbf{A}}_{\widetilde{p}}
3105 // =
3106 // \mathsf{\mathbf{F}}_{\widetilde{p}}
3107 // -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
3108 // $
3109 A.block(p_dof) += system_rhs.block(p_dof);
3110 // $
3111 // d\mathsf{\mathbf{\widetilde{J}}}
3112 // =
3113 // \mathsf{\mathbf{K}}^{-1}_{\widetilde{p}\widetilde{J}}
3114 // [
3115 // \mathsf{\mathbf{F}}_{\widetilde{p}}
3116 // -\mathsf{\mathbf{K}}_{\widetilde{p}u} d \mathsf{\mathbf{u}}
3117 // ]
3118 // $
3119 tangent_matrix.block(p_dof, J_dof)
3120 .vmult(newton_update.block(J_dof), A.block(p_dof));
3121 }
3122
3123 // we ensure here that any Dirichlet constraints
3124 // are distributed on the updated solution:
3125 constraints.distribute(newton_update);
3126
3127 // Finally we solve for the pressure
3128 // update with the substitution:
3129 // $
3130 // d \widetilde{\mathsf{\mathbf{p}}}
3131 // =
3132 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
3133 // \bigl[
3134 // \mathsf{\mathbf{F}}_{\widetilde{J}}
3135 // - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
3136 // d \widetilde{\mathsf{\mathbf{J}}}
3137 // \bigr]
3138 // $
3139 {
3140 // $
3141 // \mathsf{\mathbf{A}}_{\widetilde{J}}
3142 // =
3143 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
3144 // d \widetilde{\mathsf{\mathbf{J}}}
3145 // $
3146 tangent_matrix.block(J_dof, J_dof)
3147 .vmult(A.block(J_dof), newton_update.block(J_dof));
3148 // $
3149 // \mathsf{\mathbf{A}}_{\widetilde{J}}
3150 // =
3151 // -\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
3152 // d \widetilde{\mathsf{\mathbf{J}}}
3153 // $
3154 A.block(J_dof) *= -1.0;
3155 // $
3156 // \mathsf{\mathbf{A}}_{\widetilde{J}}
3157 // =
3158 // \mathsf{\mathbf{F}}_{\widetilde{J}}
3159 // -
3160 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
3161 // d \widetilde{\mathsf{\mathbf{J}}}
3162 // $
3163 A.block(J_dof) += system_rhs.block(J_dof);
3164 // and finally....
3165 // $
3166 // d \widetilde{\mathsf{\mathbf{p}}}
3167 // =
3168 // \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}
3169 // \bigl[
3170 // \mathsf{\mathbf{F}}_{\widetilde{J}}
3171 // - \mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{J}}
3172 // d \widetilde{\mathsf{\mathbf{J}}}
3173 // \bigr]
3174 // $
3175 tangent_matrix.block(p_dof, J_dof)
3176 .Tvmult(newton_update.block(p_dof), A.block(J_dof));
3177 }
3178
3179 // We are now at the end, so we distribute all
3180 // constrained dofs back to the Newton
3181 // update:
3182 constraints.distribute(newton_update);
3183
3184 timer.leave_subsection();
3185 }
3186 else
3187 {
3188 std::cout << " ------ " << std::flush;
3189
3190 timer.enter_subsection("Linear solver");
3191 std::cout << " SLV " << std::flush;
3192
3193 if (parameters.type_lin == "CG")
3194 {
3195 // Manual condensation of the dilatation and pressure fields on
3196 // a local level, and subsequent post-processing, took quite a
3197 // bit of effort to achieve. To recap, we had to produce the
3198 // inverse matrix
3199 // $\mathsf{\mathbf{K}}_{\widetilde{p}\widetilde{J}}^{-1}$, which
3200 // was permanently written into the global tangent matrix. We then
3201 // permanently modified $\mathsf{\mathbf{K}}_{uu}$ to produce
3202 // $\mathsf{\mathbf{K}}_{\textrm{con}}$. This involved the
3203 // extraction and manipulation of local sub-blocks of the tangent
3204 // matrix. After solving for the displacement, the individual
3205 // matrix-vector operations required to solve for dilatation and
3206 // pressure were carefully implemented. Contrast these many sequence
3207 // of steps to the much simpler and transparent implementation using
3208 // functionality provided by the LinearOperator class.
3209
3210 // For ease of later use, we define some aliases for
3211 // blocks in the RHS vector
3212 const Vector<double> &f_u = system_rhs.block(u_dof);
3213 const Vector<double> &f_p = system_rhs.block(p_dof);
3214 const Vector<double> &f_J = system_rhs.block(J_dof);
3215
3216 // ... and for blocks in the Newton update vector.
3217 Vector<double> &d_u = newton_update.block(u_dof);
3218 Vector<double> &d_p = newton_update.block(p_dof);
3219 Vector<double> &d_J = newton_update.block(J_dof);
3220
3221 // We next define some linear operators for the tangent matrix
3222 // sub-blocks We will exploit the symmetry of the system, so not all
3223 // blocks are required.
3224 const auto K_uu =
3225 linear_operator(tangent_matrix.block(u_dof, u_dof));
3226 const auto K_up =
3227 linear_operator(tangent_matrix.block(u_dof, p_dof));
3228 const auto K_pu =
3229 linear_operator(tangent_matrix.block(p_dof, u_dof));
3230 const auto K_Jp =
3231 linear_operator(tangent_matrix.block(J_dof, p_dof));
3232 const auto K_JJ =
3233 linear_operator(tangent_matrix.block(J_dof, J_dof));
3234
3235 // We then construct a LinearOperator that represents the inverse of
3236 // (square block)
3237 // $\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}$. Since it is
3238 // diagonal (or, when a higher order ansatz it used, nearly
3239 // diagonal), a Jacobi preconditioner is suitable.
3240 PreconditionSelector<SparseMatrix<double>, Vector<double>>
3241 preconditioner_K_Jp_inv("jacobi");
3242 preconditioner_K_Jp_inv.use_matrix(
3243 tangent_matrix.block(J_dof, p_dof));
3244 ReductionControl solver_control_K_Jp_inv(
3245 static_cast<unsigned int>(tangent_matrix.block(J_dof, p_dof).m() *
3246 parameters.max_iterations_lin),
3247 1.0e-30,
3248 parameters.tol_lin);
3249 SolverSelector<Vector<double>> solver_K_Jp_inv;
3250 solver_K_Jp_inv.select("cg");
3251 solver_K_Jp_inv.set_control(solver_control_K_Jp_inv);
3252 const auto K_Jp_inv =
3253 inverse_operator(K_Jp, solver_K_Jp_inv, preconditioner_K_Jp_inv);
3254
3255 // Now we can construct that transpose of
3256 // $\mathsf{\mathbf{K}}_{\widetilde{J}\widetilde{p}}^{-1}$ and a
3257 // linear operator that represents the condensed operations
3258 // $\overline{\mathsf{\mathbf{K}}}$ and
3259 // $\overline{\overline{\mathsf{\mathbf{K}}}}$ and the final
3260 // augmented matrix
3261 // $\mathsf{\mathbf{K}}_{\textrm{con}}$.
3262 // Note that the schur_complement() operator could also be of use
3263 // here, but for clarity and the purpose of demonstrating the
3264 // similarities between the formulation and implementation of the
3265 // linear solution scheme, we will perform these operations
3266 // manually.
3267 const auto K_pJ_inv = transpose_operator(K_Jp_inv);
3268 const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv;
3269 const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu;
3270 const auto K_uu_con = K_uu + K_uu_bar_bar;
3271
3272 // Lastly, we define an operator for inverse of augmented stiffness
3273 // matrix, namely $\mathsf{\mathbf{K}}_{\textrm{con}}^{-1}$. Note
3274 // that the preconditioner for the augmented stiffness matrix is
3275 // different to the case when we use static condensation. In this
3276 // instance, the preconditioner is based on a non-modified
3277 // $\mathsf{\mathbf{K}}_{uu}$, while with the first approach we
3278 // actually modified the entries of this sub-block. However, since
3279 // $\mathsf{\mathbf{K}}_{\textrm{con}}$ and
3280 // $\mathsf{\mathbf{K}}_{uu}$ operate on the same space, it remains
3281 // adequate for this problem.
3282 PreconditionSelector<SparseMatrix<double>, Vector<double>>
3283 preconditioner_K_con_inv(parameters.preconditioner_type,
3284 parameters.preconditioner_relaxation);
3285 preconditioner_K_con_inv.use_matrix(
3286 tangent_matrix.block(u_dof, u_dof));
3287 ReductionControl solver_control_K_con_inv(
3288 static_cast<unsigned int>(tangent_matrix.block(u_dof, u_dof).m() *
3289 parameters.max_iterations_lin),
3290 1.0e-30,
3291 parameters.tol_lin);
3292 SolverSelector<Vector<double>> solver_K_con_inv;
3293 solver_K_con_inv.select("cg");
3294 solver_K_con_inv.set_control(solver_control_K_con_inv);
3295 const auto K_uu_con_inv =
3296 inverse_operator(K_uu_con,
3297 solver_K_con_inv,
3298 preconditioner_K_con_inv);
3299
3300 // Now we are in a position to solve for the displacement field.
3301 // We can nest the linear operations, and the result is immediately
3302 // written to the Newton update vector.
3303 // It is clear that the implementation closely mimics the derivation
3304 // stated in the introduction.
3305 d_u =
3306 K_uu_con_inv * (f_u - K_up * (K_Jp_inv * f_J - K_pp_bar * f_p));
3307
3308 timer.leave_subsection();
3309
3310 // The operations need to post-process for the dilatation and
3311 // pressure fields are just as easy to express.
3312 timer.enter_subsection("Linear solver postprocessing");
3313 std::cout << " PP " << std::flush;
3314
3315 d_J = K_pJ_inv * (f_p - K_pu * d_u);
3316 d_p = K_Jp_inv * (f_J - K_JJ * d_J);
3317
3318 lin_it = solver_control_K_con_inv.last_step();
3319 lin_res = solver_control_K_con_inv.last_value();
3320 }
3321 else if (parameters.type_lin == "Direct")
3322 {
3323 // Solve the full block system with
3324 // a direct solver. As it is relatively
3325 // robust, it may be immune to problem
3326 // arising from the presence of the zero
3327 // $\mathsf{\mathbf{K}}_{ \widetilde{J} \widetilde{J}}$
3328 // block.
3329 SparseDirectUMFPACK A_direct;
3330 A_direct.initialize(tangent_matrix);
3331 A_direct.vmult(newton_update, system_rhs);
3332
3333 lin_it = 1;
3334 lin_res = 0.0;
3335
3336 std::cout << " -- " << std::flush;
3337 }
3338 else
3339 Assert(false, ExcMessage("Linear solver type not implemented"));
3340
3341 timer.leave_subsection();
3342
3343 // Finally, we again ensure here that any Dirichlet
3344 // constraints are distributed on the updated solution:
3345 constraints.distribute(newton_update);
3346 }
3347
3348 return std::make_pair(lin_it, lin_res);
3349 }
3350
3351 // @sect4{Solid::output_results}
3352 // Here we present how the results are written to file to be viewed
3353 // using ParaView or VisIt. The method is similar to that shown in previous
3354 // tutorials so will not be discussed in detail.
3355 template <int dim>
output_results() const3356 void Solid<dim>::output_results() const
3357 {
3358 DataOut<dim> data_out;
3359 std::vector<DataComponentInterpretation::DataComponentInterpretation>
3360 data_component_interpretation(
3361 dim, DataComponentInterpretation::component_is_part_of_vector);
3362 data_component_interpretation.push_back(
3363 DataComponentInterpretation::component_is_scalar);
3364 data_component_interpretation.push_back(
3365 DataComponentInterpretation::component_is_scalar);
3366
3367 std::vector<std::string> solution_name(dim, "displacement");
3368 solution_name.emplace_back("pressure");
3369 solution_name.emplace_back("dilatation");
3370
3371 DataOutBase::VtkFlags output_flags;
3372 output_flags.write_higher_order_cells = true;
3373 data_out.set_flags(output_flags);
3374
3375 data_out.attach_dof_handler(dof_handler);
3376 data_out.add_data_vector(solution_n,
3377 solution_name,
3378 DataOut<dim>::type_dof_data,
3379 data_component_interpretation);
3380
3381 // Since we are dealing with a large deformation problem, it would be nice
3382 // to display the result on a displaced grid! The MappingQEulerian class
3383 // linked with the DataOut class provides an interface through which this
3384 // can be achieved without physically moving the grid points in the
3385 // Triangulation object ourselves. We first need to copy the solution to
3386 // a temporary vector and then create the Eulerian mapping. We also
3387 // specify the polynomial degree to the DataOut object in order to produce
3388 // a more refined output data set when higher order polynomials are used.
3389 Vector<double> soln(solution_n.size());
3390 for (unsigned int i = 0; i < soln.size(); ++i)
3391 soln(i) = solution_n(i);
3392 MappingQEulerian<dim> q_mapping(degree, dof_handler, soln);
3393 data_out.build_patches(q_mapping, degree);
3394
3395 std::ofstream output("solution-" + std::to_string(dim) + "d-" +
3396 std::to_string(time.get_timestep()) + ".vtu");
3397 data_out.write_vtu(output);
3398 }
3399
3400 } // namespace Step44
3401
3402
3403 // @sect3{Main function}
3404 // Lastly we provide the main driver function which appears
3405 // no different to the other tutorials.
main()3406 int main()
3407 {
3408 using namespace Step44;
3409
3410 try
3411 {
3412 const unsigned int dim = 3;
3413 Solid<dim> solid("parameters.prm");
3414 solid.run();
3415 }
3416 catch (std::exception &exc)
3417 {
3418 std::cerr << std::endl
3419 << std::endl
3420 << "----------------------------------------------------"
3421 << std::endl;
3422 std::cerr << "Exception on processing: " << std::endl
3423 << exc.what() << std::endl
3424 << "Aborting!" << std::endl
3425 << "----------------------------------------------------"
3426 << std::endl;
3427
3428 return 1;
3429 }
3430 catch (...)
3431 {
3432 std::cerr << std::endl
3433 << std::endl
3434 << "----------------------------------------------------"
3435 << std::endl;
3436 std::cerr << "Unknown exception!" << std::endl
3437 << "Aborting!" << std::endl
3438 << "----------------------------------------------------"
3439 << std::endl;
3440 return 1;
3441 }
3442
3443 return 0;
3444 }
3445