From 821feaf203604b5038f53206f2d8c425e4636a3d Mon Sep 17 00:00:00 2001 From: Jean-Paul Pelteret Date: Mon, 21 Dec 2015 10:41:46 +0100 Subject: [PATCH] A classical approach to solving the cook membrane problem. In this work we take a classical approach to solving the equations governing quasi-static finite-strain compressible elasticity, with code based on `step-44`. The formulation adopted here is that seen in many texts on solid mechanics and can be used as the starting point for extension into many topics such as material anisotropy, rate dependence or plasticity, or even as a component of multi-physics problems. --- .../CMakeLists.txt | 36 + .../README.md | 117 + .../cook_membrane.cc | 2370 +++++++++++++++++ .../doc/author | 2 + .../doc/builds-on | 1 + .../doc/dependencies | 0 .../doc/disp_soln.png | Bin 0 -> 28134 bytes .../doc/problem_setup.png | Bin 0 -> 17787 bytes .../doc/ref_grid.png | Bin 0 -> 14414 bytes .../doc/tooltip | 1 + 10 files changed, 2527 insertions(+) create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/README.md create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/author create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/dependencies create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/disp_soln.png create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/problem_setup.png create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/ref_grid.png create mode 100644 Quasi_static_Finite_strain_Compressible_Elasticity/doc/tooltip diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt new file mode 100644 index 0000000..724337c --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/CMakeLists.txt @@ -0,0 +1,36 @@ +## +# CMake script for the step- tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "cook_membrane") + +# Declare all source files the target consists of: +SET(TARGET_SRC + ${TARGET}.cc + # You can specify additional files here! + ) + +SET(CLEAN_UP_FILES + # a custom list of globs, e.g. *.log *.vtk + *.vtk +) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.0 QUIET + HINTS ${deal.II_DIR} ${DEAL_II_DIR} ../ ../../ $ENV{DEAL_II_DIR} + ) +IF(NOT ${deal.II_FOUND}) + MESSAGE(FATAL_ERROR "\n" + "*** Could not locate deal.II. ***\n\n" + "You may want to either pass a flag -DDEAL_II_DIR=/path/to/deal.II to cmake\n" + "or set an environment variable \"DEAL_II_DIR\" that contains this path." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/README.md b/Quasi_static_Finite_strain_Compressible_Elasticity/README.md new file mode 100644 index 0000000..3605f01 --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/README.md @@ -0,0 +1,117 @@ +# Gallery Entry: deal.II-Quasi-static_Finite_strain_Compressible_Elasticity + + +## Overview +The Cook membrane (or cantilever) problem is a classic benchmark test for +finite element formulations for solid mechanics. It is typically used to +test for and demonstrate the shear-locking (or locking-free) behaviour of +a finite element ansatz under quasi-incompressible conditions. +As it is so widely referred to in the literature on finite-strain elasticity, +we reproduce the example here. +However, we consider on the compressible case to avoid many of the complexities +that arise in `step-44`, which provides an efficient approach to deal with +the quasi-incompressible case. + +### A classical approach to solving the cook membrane problem. + +In this work we take a classical approach to solving the equations governing +quasi-static finite-strain compressible elasticity, with code based on +`step-44`. The formulation adopted here is that seen in many texts on solid +mechanics and can be used as the starting point for extension into many +topics such as material anisotropy, rate dependence or plasticity, or even as +a component of multi-physics problems. + +The basic problem configuration is summarised in the following image. +A beam of specific dimensions is fixed at one end and a uniform traction load +is applied at the other end such that the total force acting on this surface +totals 1 Newton. Displacement in the third coordinate direction (out of plane) +is prevented in order to impose plane strain conditions. + +![Problem geometry](./doc/problem_setup.png) + +Note that we perform a three-dimensional computation as two-dimensional +elasticity corresponds to neither plane-strain nor plane-stress conditions. + + +## Requirements +* Version 8.2.1 or greater of `deal.II` + + +## Compiling and running +Similar to the example programs, run +``` +cmake -DDEAL_II_DIR=/path/to/deal.II . +``` +in this directory to configure the problem. +You can switch between debug and release mode by calling either +``` +make debug +``` +or +``` +make release +``` +The problem may then be run with +``` +make run +``` + + +## Recommended Literature +* C. Miehe (1994), Aspects of the formulation and finite element implementation of large strain isotropic elasticity International Journal for Numerical Methods in Engineering 37 , 12, 1981-2004; +* G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons; +* P. Wriggers (2008), Nonlinear finite element methods, Springer; +* T.J.R. Hughes (2000), The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover. + +The derivation of the finite-element problem, namely the definition and +linearisation of the residual and their subsequent discretisation are quite +lengthy and involved. Thankfully, the classical approach adopted in this work is +well documented and therefore does not need to be reproduced here. +We refer the reader to, among many other possible texts, Holzapfel (2001) and +Wriggers (2008) for a detailed description of the approach applied in this work. +It amounts to a reduction and slight reworking of `step-44` (accounting for +the removal of the two additional fields used therein). We also refer the reader +to `step-44` for a brief overview of the continuum mechanics and kinematics +related to solid mechanics. + +## Results +These results were produced using the following material properties: +* Shear modulus is 422.5kPa +* Poisson ratio is 0.3 + +The 32x32x1 discretised reference geometry looks as follows: + +![Problem geometry](./doc/ref_grid.png) + +And an example of the displaced solution is given in the next image. + +![Displaced solution](./doc/disp_soln.png) + +Below we briefly document the tip displacement as predicted for different +discretisation levels and ansatz for the displacement field. +A direct comparison of the following results can be made with those found +in Miehe (1994). + +#### Number of degrees of freedom + +Elements per edge | Q1 | Q2 +:---------------: | :-------------: | :-------------: +1 | 24 | 81 +2 | 54 | 225 +4 | 150 | 729 +8 | 486 | 2601 +16 | 1734 | 9801 +32 | 6534 | 38025 +64 | 25350 | 149769 + +#### Tip y-displacement (in mm) + +Elements per edge | Q1 | Q2 +:---------------: | :-------------: | :-------------: +1 | 5.15 | 12.19 +2 | 8.72 | 13.83 +4 | 12.02 | 14.22 +8 | 13.61 | 14.30 +16 | 14.13 | 14.32 +32 | 14.28 | 14.33 +64 | 14.32 | 14.33 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc new file mode 100644 index 0000000..6dffb8b --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/cook_membrane.cc @@ -0,0 +1,2370 @@ +/* --------------------------------------------------------------------- + * Copyright (C) 2010 - 2015 by the deal.II authors and + * Jean-Paul Pelteret and Andrew McBride + * + * This file is part of the deal.II library. + * + * The deal.II library is free software; you can use it, redistribute + * it, and/or modify it under the terms of the GNU Lesser General + * Public License as published by the Free Software Foundation; either + * version 2.1 of the License, or (at your option) any later version. + * The full text of the license can be found in the file LICENSE at + * the top level of the deal.II distribution. + * + * --------------------------------------------------------------------- + */ + +/* + * Authors: Jean-Paul Pelteret, University of Erlangen-Nuremberg, + * Andrew McBride, University of Cape Town, 2015 + */ + + +// We start by including all the necessary deal.II header files and some C++ +// related ones. They have been discussed in detail in previous tutorial +// programs, so you need only refer to past tutorials for details. +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include +#include +#include +#include // TEMP + +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include + + +// We then stick everything that relates to this tutorial program into a +// namespace of its own, and import all the deal.II function and class names +// into it: +namespace Cook_Membrane +{ + using namespace dealii; + +// @sect3{Run-time parameters} +// +// There are several parameters that can be set in the code so we set up a +// ParameterHandler object to read in the choices at run-time. + namespace Parameters + { +// @sect4{Finite Element system} + +// Here we specify the polynomial order used to approximate the solution. +// The quadrature order should be adjusted accordingly. + struct FESystem + { + unsigned int poly_degree; + unsigned int quad_order; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + + void FESystem::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + prm.declare_entry("Polynomial degree", "2", + Patterns::Integer(0), + "Displacement system polynomial order"); + + prm.declare_entry("Quadrature order", "3", + Patterns::Integer(0), + "Gauss quadrature order"); + } + prm.leave_subsection(); + } + + void FESystem::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Finite element system"); + { + poly_degree = prm.get_integer("Polynomial degree"); + quad_order = prm.get_integer("Quadrature order"); + } + prm.leave_subsection(); + } + +// @sect4{Geometry} + +// Make adjustments to the problem geometry and its discretisation. + struct Geometry + { + unsigned int elements_per_edge; + double scale; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Geometry::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + prm.declare_entry("Elements per edge", "32", + Patterns::Integer(0), + "Number of elements per long edge of the beam"); + + prm.declare_entry("Grid scale", "1e-3", + Patterns::Double(0.0), + "Global grid scaling factor"); + } + prm.leave_subsection(); + } + + void Geometry::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + elements_per_edge = prm.get_integer("Elements per edge"); + scale = prm.get_double("Grid scale"); + } + prm.leave_subsection(); + } + +// @sect4{Materials} + +// We also need the shear modulus $ \mu $ and Poisson ration $ \nu $ for the +// neo-Hookean material. + struct Materials + { + double nu; + double mu; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Materials::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + prm.declare_entry("Poisson's ratio", "0.4999", + Patterns::Double(-1.0,0.5), + "Poisson's ratio"); + + prm.declare_entry("Shear modulus", "80.194e6", + Patterns::Double(), + "Shear modulus"); + } + prm.leave_subsection(); + } + + void Materials::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + nu = prm.get_double("Poisson's ratio"); + mu = prm.get_double("Shear modulus"); + } + prm.leave_subsection(); + } + +// @sect4{Linear solver} + +// Next, we choose both solver and preconditioner settings. The use of an +// effective preconditioner is critical to ensure convergence when a large +// nonlinear motion occurs within a Newton increment. + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + std::string preconditioner_type; + double preconditioner_relaxation; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void LinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + prm.declare_entry("Solver type", "CG", + Patterns::Selection("CG|Direct"), + "Type of solver used to solve the linear system"); + + prm.declare_entry("Residual", "1e-6", + Patterns::Double(0.0), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry("Max iteration multiplier", "1", + Patterns::Double(0.0), + "Linear solver iterations (multiples of the system matrix size)"); + + prm.declare_entry("Preconditioner type", "ssor", + Patterns::Selection("jacobi|ssor"), + "Type of preconditioner"); + + prm.declare_entry("Preconditioner relaxation", "0.65", + Patterns::Double(0.0), + "Preconditioner relaxation value"); + } + prm.leave_subsection(); + } + + void LinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Linear solver"); + { + type_lin = prm.get("Solver type"); + tol_lin = prm.get_double("Residual"); + max_iterations_lin = prm.get_double("Max iteration multiplier"); + preconditioner_type = prm.get("Preconditioner type"); + preconditioner_relaxation = prm.get_double("Preconditioner relaxation"); + } + prm.leave_subsection(); + } + +// @sect4{Nonlinear solver} + +// A Newton-Raphson scheme is used to solve the nonlinear system of governing +// equations. We now define the tolerances and the maximum number of +// iterations for the Newton-Raphson nonlinear solver. + struct NonlinearSolver + { + unsigned int max_iterations_NR; + double tol_f; + double tol_u; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void NonlinearSolver::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + prm.declare_entry("Max iterations Newton-Raphson", "10", + Patterns::Integer(0), + "Number of Newton-Raphson iterations allowed"); + + prm.declare_entry("Tolerance force", "1.0e-9", + Patterns::Double(0.0), + "Force residual tolerance"); + + prm.declare_entry("Tolerance displacement", "1.0e-6", + Patterns::Double(0.0), + "Displacement error tolerance"); + } + prm.leave_subsection(); + } + + void NonlinearSolver::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Nonlinear solver"); + { + max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson"); + tol_f = prm.get_double("Tolerance force"); + tol_u = prm.get_double("Tolerance displacement"); + } + prm.leave_subsection(); + } + +// @sect4{Time} + +// Set the timestep size $ \varDelta t $ and the simulation end-time. + struct Time + { + double delta_t; + double end_time; + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + void Time::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + prm.declare_entry("End time", "1", + Patterns::Double(), + "End time"); + + prm.declare_entry("Time step size", "0.1", + Patterns::Double(), + "Time step size"); + } + prm.leave_subsection(); + } + + void Time::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Time"); + { + end_time = prm.get_double("End time"); + delta_t = prm.get_double("Time step size"); + } + prm.leave_subsection(); + } + +// @sect4{All parameters} + +// Finally we consolidate all of the above structures into a single container +// that holds all of our run-time selections. + struct AllParameters : public FESystem, + public Geometry, + public Materials, + public LinearSolver, + public NonlinearSolver, + public Time + + { + AllParameters(const std::string &input_file); + + static void + declare_parameters(ParameterHandler &prm); + + void + parse_parameters(ParameterHandler &prm); + }; + + AllParameters::AllParameters(const std::string &input_file) + { + ParameterHandler prm; + declare_parameters(prm); + prm.read_input(input_file); + parse_parameters(prm); + } + + void AllParameters::declare_parameters(ParameterHandler &prm) + { + FESystem::declare_parameters(prm); + Geometry::declare_parameters(prm); + Materials::declare_parameters(prm); + LinearSolver::declare_parameters(prm); + NonlinearSolver::declare_parameters(prm); + Time::declare_parameters(prm); + } + + void AllParameters::parse_parameters(ParameterHandler &prm) + { + FESystem::parse_parameters(prm); + Geometry::parse_parameters(prm); + Materials::parse_parameters(prm); + LinearSolver::parse_parameters(prm); + NonlinearSolver::parse_parameters(prm); + Time::parse_parameters(prm); + } + } + +// @sect3{Some standard tensors} + +// Now we define some frequently used second and fourth-order tensors: + template + class StandardTensors + { + public: + + // $\mathbf{I}$ + static const SymmetricTensor<2, dim> I; + // $\mathbf{I} \otimes \mathbf{I}$ + static const SymmetricTensor<4, dim> IxI; + // $\mathcal{S}$, note that as we only use this fourth-order unit tensor + // to operate on symmetric second-order tensors. To maintain notation + // consistent with Holzapfel (2001) we name the tensor $\mathcal{I}$ + static const SymmetricTensor<4, dim> II; + // Fourth-order deviatoric tensor such that + // $\textrm{dev} \{ \bullet \} = \{ \bullet \} - + // [1/\textrm{dim}][ \{ \bullet\} :\mathbf{I}]\mathbf{I}$ + static const SymmetricTensor<4, dim> dev_P; + }; + + template + const SymmetricTensor<2, dim> + StandardTensors::I = unit_symmetric_tensor(); + + template + const SymmetricTensor<4, dim> + StandardTensors::IxI = outer_product(I, I); + + template + const SymmetricTensor<4, dim> + StandardTensors::II = identity_tensor(); + + template + const SymmetricTensor<4, dim> + StandardTensors::dev_P = deviator_tensor(); + +// @sect3{Time class} + +// A simple class to store time data. Its functioning is transparent so no +// discussion is necessary. For simplicity we assume a constant time step +// size. + class Time + { + public: + Time (const double time_end, + const double delta_t) + : + timestep(0), + time_current(0.0), + time_end(time_end), + delta_t(delta_t) + {} + + virtual ~Time() + {} + + double current() const + { + return time_current; + } + double end() const + { + return time_end; + } + double get_delta_t() const + { + return delta_t; + } + unsigned int get_timestep() const + { + return timestep; + } + void increment() + { + time_current += delta_t; + ++timestep; + } + + private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; + }; + +// @sect3{Compressible neo-Hookean material within a one-field formulation} + +// As discussed in the literature and step-44, Neo-Hookean materials are a type +// of hyperelastic materials. The entire domain is assumed to be composed of a +// compressible neo-Hookean material. This class defines the behaviour of +// this material within a one-field formulation. Compressible neo-Hookean +// materials can be described by a strain-energy function (SEF) $ \Psi = +// \Psi_{\text{iso}}(\overline{\mathbf{b}}) + \Psi_{\text{vol}}(J) +// $. +// +// The isochoric response is given by $ +// \Psi_{\text{iso}}(\overline{\mathbf{b}}) = c_{1} [\overline{I}_{1} - 3] $ +// where $ c_{1} = \frac{\mu}{2} $ and $\overline{I}_{1}$ is the first +// invariant of the left- or right-isochoric Cauchy-Green deformation tensors. +// That is $\overline{I}_1 :=\textrm{tr}(\overline{\mathbf{b}})$. In this +// example the SEF that governs the volumetric response is defined as $ +// \Psi_{\text{vol}}(J) = \kappa \frac{1}{4} [ J^2 - 1 +// - 2\textrm{ln}\; J ]$, where $\kappa:= \lambda + 2/3 \mu$ is +// the bulk modulus +// and $\lambda$ is Lame's first +// parameter. +// +// The following class will be used to characterize the material we work with, +// and provides a central point that one would need to modify if one were to +// implement a different material model. For it to work, we will store one +// object of this type per quadrature point, and in each of these objects +// store the current state (characterized by the values or measures of the +// displacement field) so that we can compute the elastic coefficients +// linearized around the current state. + template + class Material_Compressible_Neo_Hook_One_Field + { + public: + Material_Compressible_Neo_Hook_One_Field(const double mu, + const double nu) + : + kappa((2.0 * mu * (1.0 + nu)) / (3.0 * (1.0 - 2.0 * nu))), + c_1(mu / 2.0), + det_F(1.0), + b_bar(StandardTensors::I) + { + Assert(kappa > 0, ExcInternalError()); + } + + ~Material_Compressible_Neo_Hook_One_Field() + {} + + // We update the material model with various deformation dependent data + // based on the deformation gradient $\mathbf{F}$ and the volumetric Jacobian + // $J = \text{det} \mathbf{F}$, and at the end of the function include a + // physical check for internal consistency: + void update_material_data(const Tensor<2, dim> &F) + { + det_F = determinant(F); + b_bar = std::pow(det_F, -2.0 / 3.0) * symmetrize(F * transpose(F)); + + Assert(det_F > 0, ExcInternalError()); + } + + // The second function determines the Kirchhoff stress $\boldsymbol{\tau} + // = \boldsymbol{\tau}_{\textrm{iso}} + \boldsymbol{\tau}_{\textrm{vol}}$ + SymmetricTensor<2, dim> get_tau() + { + // See Holzapfel p231 eq6.98 onwards + return get_tau_iso() + get_tau_vol(); + } + + // The fourth-order elasticity tensor in the spatial setting + // $\mathfrak{c}$ is calculated from the SEF $\Psi$ as $ J + // \mathfrak{c}_{ijkl} = F_{iA} F_{jB} \mathfrak{C}_{ABCD} F_{kC} F_{lD}$ + // where $ \mathfrak{C} = 4 \frac{\partial^2 \Psi(\mathbf{C})}{\partial + // \mathbf{C} \partial \mathbf{C}}$ + SymmetricTensor<4, dim> get_Jc() const + { + return get_Jc_vol() + get_Jc_iso(); + } + + // Derivative of the volumetric free energy with respect to + // $J$ return $\frac{\partial + // \Psi_{\text{vol}}(J)}{\partial J}$ + double get_dPsi_vol_dJ() const + { + return (kappa / 2.0) * (det_F - 1.0 / det_F); + } + + // Second derivative of the volumetric free energy wrt $J$. We + // need the following computation explicitly in the tangent so we make it + // public. We calculate $\frac{\partial^2 + // \Psi_{\textrm{vol}}(J)}{\partial J \partial + // J}$ + double get_d2Psi_vol_dJ2() const + { + return ( (kappa / 2.0) * (1.0 + 1.0 / (det_F * det_F))); + } + + // The next few functions return various data that we choose to store with + // the material: + double get_det_F() const + { + return det_F; + } + + protected: + // Define constitutive model parameters $\kappa$ (bulk modulus) and the + // neo-Hookean model parameter $c_1$: + const double kappa; + const double c_1; + + // Model specific data that is convenient to store with the material: + double det_F; + SymmetricTensor<2, dim> b_bar; + + // The following functions are used internally in determining the result + // of some of the public functions above. The first one determines the + // volumetric Kirchhoff stress $\boldsymbol{\tau}_{\textrm{vol}}$. + // Note the difference in its definition when compared to step-44. + SymmetricTensor<2, dim> get_tau_vol() const + { + return get_dPsi_vol_dJ() * det_F * StandardTensors::I; + } + + // Next, determine the isochoric Kirchhoff stress + // $\boldsymbol{\tau}_{\textrm{iso}} = + // \mathcal{P}:\overline{\boldsymbol{\tau}}$: + SymmetricTensor<2, dim> get_tau_iso() const + { + return StandardTensors::dev_P * get_tau_bar(); + } + + // Then, determine the fictitious Kirchhoff stress + // $\overline{\boldsymbol{\tau}}$: + SymmetricTensor<2, dim> get_tau_bar() const + { + return 2.0 * c_1 * b_bar; + } + + // Calculate the volumetric part of the tangent $J + // \mathfrak{c}_\textrm{vol}$. Again, note the difference in its + // definition when compared to step-44. The extra terms result from two + // quantities in $\boldsymbol{\tau}_{\textrm{vol}}$ being dependent on + // $\boldsymbol{F}$. + SymmetricTensor<4, dim> get_Jc_vol() const + { + // See Holzapfel p265 + return det_F + * ( (get_dPsi_vol_dJ() + det_F * get_d2Psi_vol_dJ2())*StandardTensors::IxI + - (2.0 * get_dPsi_vol_dJ())*StandardTensors::II ); + } + + // Calculate the isochoric part of the tangent $J + // \mathfrak{c}_\textrm{iso}$: + SymmetricTensor<4, dim> get_Jc_iso() const + { + const SymmetricTensor<2, dim> tau_bar = get_tau_bar(); + const SymmetricTensor<2, dim> tau_iso = get_tau_iso(); + const SymmetricTensor<4, dim> tau_iso_x_I + = outer_product(tau_iso, + StandardTensors::I); + const SymmetricTensor<4, dim> I_x_tau_iso + = outer_product(StandardTensors::I, + tau_iso); + const SymmetricTensor<4, dim> c_bar = get_c_bar(); + + return (2.0 / 3.0) * trace(tau_bar) + * StandardTensors::dev_P + - (2.0 / 3.0) * (tau_iso_x_I + I_x_tau_iso) + + StandardTensors::dev_P * c_bar + * StandardTensors::dev_P; + } + + // Calculate the fictitious elasticity tensor $\overline{\mathfrak{c}}$. + // For the material model chosen this is simply zero: + SymmetricTensor<4, dim> get_c_bar() const + { + return SymmetricTensor<4, dim>(); + } + }; + +// @sect3{Quadrature point history} + +// As seen in step-18, the PointHistory class offers a method +// for storing data at the quadrature points. Here each quadrature point +// holds a pointer to a material description. Thus, different material models +// can be used in different regions of the domain. Among other data, we +// choose to store the Kirchhoff stress $\boldsymbol{\tau}$ and the tangent +// $J\mathfrak{c}$ for the quadrature points. + template + class PointHistory + { + public: + PointHistory() + : + material(NULL), + F_inv(StandardTensors::I), + tau(SymmetricTensor<2, dim>()), + Jc(SymmetricTensor<4, dim>()) + {} + + virtual ~PointHistory() + { + delete material; + material = NULL; + } + + // The first function is used to create a material object and to + // initialize all tensors correctly: The second one updates the stored + // values and stresses based on the current deformation measure + // $\textrm{Grad}\mathbf{u}_{\textrm{n}}$. + void setup_lqp (const Parameters::AllParameters ¶meters) + { + material = new Material_Compressible_Neo_Hook_One_Field(parameters.mu, + parameters.nu); + update_values(Tensor<2, dim>()); + } + + // To this end, we calculate the deformation gradient $\mathbf{F}$ from + // the displacement gradient $\textrm{Grad}\ \mathbf{u}$, i.e. + // $\mathbf{F}(\mathbf{u}) = \mathbf{I} + \textrm{Grad}\ \mathbf{u}$ and + // then let the material model associated with this quadrature point + // update itself. When computing the deformation gradient, we have to take + // care with which data types we compare the sum $\mathbf{I} + + // \textrm{Grad}\ \mathbf{u}$: Since $I$ has data type SymmetricTensor, + // just writing I + Grad_u_n would convert the second + // argument to a symmetric tensor, perform the sum, and then cast the + // result to a Tensor (i.e., the type of a possibly nonsymmetric + // tensor). However, since Grad_u_n is nonsymmetric in + // general, the conversion to SymmetricTensor will fail. We can avoid this + // back and forth by converting $I$ to Tensor first, and then performing + // the addition as between nonsymmetric tensors: + void update_values (const Tensor<2, dim> &Grad_u_n) + { + const Tensor<2, dim> F + = (Tensor<2, dim>(StandardTensors::I) + + Grad_u_n); + material->update_material_data(F); + + // The material has been updated so we now calculate the Kirchhoff + // stress $\mathbf{\tau}$, the tangent $J\mathfrak{c}$ and the first and + // second derivatives of the volumetric free energy. + // + // We also store the inverse of the deformation gradient since we + // frequently use it: + F_inv = invert(F); + tau = material->get_tau(); + Jc = material->get_Jc(); + } + + // We offer an interface to retrieve certain data. Here are the kinematic + // variables: + + double get_det_F() const + { + return material->get_det_F(); + } + + const Tensor<2, dim> &get_F_inv() const + { + return F_inv; + } + + // ...and the kinetic variables. These are used in the material and + // global tangent matrix and residual assembly operations: + + const SymmetricTensor<2, dim> &get_tau() const + { + return tau; + } + + // And finally the tangent: + const SymmetricTensor<4, dim> &get_Jc() const + { + return Jc; + } + + // In terms of member functions, this class stores for the quadrature + // point it represents a copy of a material type in case different + // materials are used in different regions of the domain, as well as the + // inverse of the deformation gradient... + private: + Material_Compressible_Neo_Hook_One_Field *material; + + Tensor<2, dim> F_inv; + + // ... and stress-type variables along with the tangent $J\mathfrak{c}$: + SymmetricTensor<2, dim> tau; + SymmetricTensor<4, dim> Jc; + }; + + +// @sect3{Quasi-static compressible finite-strain solid} + +// The Solid class is the central class in that it represents the problem at +// hand. It follows the usual scheme in that all it really has is a +// constructor, destructor and a run() function that dispatches +// all the work to private functions of this class: + template + class Solid + { + public: + Solid(const std::string &input_file); + + virtual + ~Solid(); + + void + run(); + + private: + + // In the private section of this class, we first forward declare a number + // of objects that are used in parallelizing work using the WorkStream + // object (see the @ref threads module for more information on this). + // + // We declare such structures for the computation of tangent (stiffness) + // matrix, right hand side, and for updating quadrature points: + struct PerTaskData_K; + struct ScratchData_K; + + struct PerTaskData_RHS; + struct ScratchData_RHS; + + struct PerTaskData_UQPH; + struct ScratchData_UQPH; + + // We start the collection of member functions with one that builds the + // grid: + void + make_grid(); + + // Set up the finite element system to be solved: + void + system_setup(); + + // Several functions to assemble the system and right hand side matrices + // using multithreading. Each of them comes as a wrapper function, one + // that is executed to do the work in the WorkStream model on one cell, + // and one that copies the work done on this one cell into the global + // object that represents it: + void + assemble_system_tangent(); + + void + assemble_system_tangent_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K &scratch, + PerTaskData_K &data); + + void + copy_local_to_global_K(const PerTaskData_K &data); + + void + assemble_system_rhs(); + + void + assemble_system_rhs_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS &scratch, + PerTaskData_RHS &data); + + void + copy_local_to_global_rhs(const PerTaskData_RHS &data); + + // Apply Dirichlet boundary conditions on the displacement field + void + make_constraints(const int &it_nr); + + // Create and update the quadrature points. Here, no data needs to be + // copied into a global object, so the copy_local_to_global function is + // empty: + void + setup_qph(); + + void + update_qph_incremental(const BlockVector &solution_delta); + + void + update_qph_incremental_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH &scratch, + PerTaskData_UQPH &data); + + void + copy_local_to_global_UQPH(const PerTaskData_UQPH &/*data*/) + {} + + // Solve for the displacement using a Newton-Raphson method. We break this + // function into the nonlinear loop and the function that solves the + // linearized Newton-Raphson step: + void + solve_nonlinear_timestep(BlockVector &solution_delta); + + std::pair + solve_linear_system(BlockVector &newton_update); + + // Solution retrieval as well as post-processing and writing data to file: + BlockVector + get_total_solution(const BlockVector &solution_delta) const; + + void + output_results() const; + + // Finally, some member variables that describe the current state: A + // collection of the parameters used to describe the problem setup... + Parameters::AllParameters parameters; + + // ...the volume of the reference and current configurations... + double vol_reference; + double vol_current; + + // ...and description of the geometry on which the problem is solved: + Triangulation triangulation; + + // Also, keep track of the current time and the time spent evaluating + // certain functions + Time time; + TimerOutput timer; + + // A storage object for quadrature point information. See step-18 for + // more on this: + std::vector > quadrature_point_history; + + // A description of the finite-element system including the displacement + // polynomial degree, the degree-of-freedom handler, number of DoFs per + // cell and the extractor objects used to retrieve information from the + // solution vectors: + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler_ref; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + + // Description of how the block-system is arranged. There is just 1 block, + // that contains a vector DOF $\mathbf{u}$. + // There are two reasons that we retain the block system in this problem. + // The first is pure laziness to perform further modifications to the + // code from which this work originated. The second is that a block system + // would typically necessary when extending this code to multiphysics + // problems. + static const unsigned int n_blocks = 1; + static const unsigned int n_components = dim; + static const unsigned int first_u_component = 0; + + enum + { + u_dof = 0 + }; + + std::vector dofs_per_block; + + // Rules for Gauss-quadrature on both the cell and faces. The number of + // quadrature points on both cells and faces is recorded. + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + + // Objects that store the converged solution and right-hand side vectors, + // as well as the tangent matrix. There is a ConstraintMatrix object used + // to keep track of constraints. We make use of a sparsity pattern + // designed for a block system. + ConstraintMatrix constraints; + BlockSparsityPattern sparsity_pattern; + BlockSparseMatrix tangent_matrix; + BlockVector system_rhs; + BlockVector solution_n; + + // Then define a number of variables to store norms and update norms and + // normalisation factors. + struct Errors + { + Errors() + : + norm(1.0), u(1.0) + {} + + void reset() + { + norm = 1.0; + u = 1.0; + } + void normalise(const Errors &rhs) + { + if (rhs.norm != 0.0) + norm /= rhs.norm; + if (rhs.u != 0.0) + u /= rhs.u; + } + + double norm, u; + }; + + Errors error_residual, error_residual_0, error_residual_norm, error_update, + error_update_0, error_update_norm; + + // Methods to calculate error measures + void + get_error_residual(Errors &error_residual); + + void + get_error_update(const BlockVector &newton_update, + Errors &error_update); + + // Print information to screen in a pleasing way... + static + void + print_conv_header(); + + void + print_conv_footer(); + + void + print_vertical_tip_displacement(); + }; + +// @sect3{Implementation of the Solid class} + +// @sect4{Public interface} + +// We initialise the Solid class using data extracted from the parameter file. + template + Solid::Solid(const std::string &input_file) + : + parameters(input_file), + triangulation(Triangulation::maximum_smoothing), + time(parameters.end_time, parameters.delta_t), + timer(std::cout, + TimerOutput::summary, + TimerOutput::wall_times), + degree(parameters.poly_degree), + // The Finite Element System is composed of dim continuous displacement + // DOFs. + fe(FE_Q(parameters.poly_degree), dim), // displacement + dof_handler_ref(triangulation), + dofs_per_cell (fe.dofs_per_cell), + u_fe(first_u_component), + dofs_per_block(n_blocks), + qf_cell(parameters.quad_order), + qf_face(parameters.quad_order), + n_q_points (qf_cell.size()), + n_q_points_f (qf_face.size()) + { + + } + +// The class destructor simply clears the data held by the DOFHandler + template + Solid::~Solid() + { + dof_handler_ref.clear(); + } + + +// In solving the quasi-static problem, the time becomes a loading parameter, +// i.e. we increasing the loading linearly with time, making the two concepts +// interchangeable. We choose to increment time linearly using a constant time +// step size. +// +// We start the function with preprocessing, and then output the initial grid +// before starting the simulation proper with the first time (and loading) +// increment. +// + template + void Solid::run() + { + make_grid(); + system_setup(); + output_results(); + time.increment(); + + // We then declare the incremental solution update $\varDelta + // \mathbf{\Xi}:= \{\varDelta \mathbf{u}\}$ and start the loop over the + // time domain. + // + // At the beginning, we reset the solution update for this time step... + BlockVector solution_delta(dofs_per_block); + while (time.current() <= time.end()) + { + solution_delta = 0.0; + + // ...solve the current time step and update total solution vector + // $\mathbf{\Xi}_{\textrm{n}} = \mathbf{\Xi}_{\textrm{n-1}} + + // \varDelta \mathbf{\Xi}$... + solve_nonlinear_timestep(solution_delta); + solution_n += solution_delta; + + // ...and plot the results before moving on happily to the next time + // step: + output_results(); + time.increment(); + } + + // Lastly, we print the vertical tip displacement of the Cook cantilever + // after the full load is applied + print_vertical_tip_displacement(); + } + + +// @sect3{Private interface} + +// @sect4{Threading-building-blocks structures} + +// The first group of private member functions is related to parallization. +// We use the Threading Building Blocks library (TBB) to perform as many +// computationally intensive distributed tasks as possible. In particular, we +// assemble the tangent matrix and right hand side vector, and update data +// stored at the quadrature points using TBB. Our main tool for this is the +// WorkStream class (see the @ref threads module for more information). + +// Firstly we deal with the tangent matrix assembly structures. The +// PerTaskData object stores local contributions. + template + struct Solid::PerTaskData_K + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + + PerTaskData_K(const unsigned int dofs_per_cell) + : + cell_matrix(dofs_per_cell, dofs_per_cell), + local_dof_indices(dofs_per_cell) + {} + + void reset() + { + cell_matrix = 0.0; + } + }; + + +// On the other hand, the ScratchData object stores the larger objects such as +// the shape-function values array (Nx) and a shape function +// gradient and symmetric gradient vector which we will use during the +// assembly. + template + struct Solid::ScratchData_K + { + FEValues fe_values_ref; + + std::vector > Nx; + std::vector > > grad_Nx; + std::vector > > symm_grad_Nx; + + ScratchData_K(const FiniteElement &fe_cell, + const QGauss &qf_cell, + const UpdateFlags uf_cell) + : + fe_values_ref(fe_cell, qf_cell, uf_cell), + Nx(qf_cell.size(), + std::vector(fe_cell.dofs_per_cell)), + grad_Nx(qf_cell.size(), + std::vector >(fe_cell.dofs_per_cell)), + symm_grad_Nx(qf_cell.size(), + std::vector > + (fe_cell.dofs_per_cell)) + {} + + ScratchData_K(const ScratchData_K &rhs) + : + fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()), + Nx(rhs.Nx), + grad_Nx(rhs.grad_Nx), + symm_grad_Nx(rhs.symm_grad_Nx) + {} + + void reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert( grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + + }; + +// Next, the same approach is used for the right-hand side assembly. The +// PerTaskData object again stores local contributions and the ScratchData +// object the shape function object and precomputed values vector: + template + struct Solid::PerTaskData_RHS + { + Vector cell_rhs; + std::vector local_dof_indices; + + PerTaskData_RHS(const unsigned int dofs_per_cell) + : + cell_rhs(dofs_per_cell), + local_dof_indices(dofs_per_cell) + {} + + void reset() + { + cell_rhs = 0.0; + } + }; + + + template + struct Solid::ScratchData_RHS + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + + std::vector > Nx; + std::vector > > symm_grad_Nx; + + ScratchData_RHS(const FiniteElement &fe_cell, + const QGauss &qf_cell, const UpdateFlags uf_cell, + const QGauss & qf_face, const UpdateFlags uf_face) + : + fe_values_ref(fe_cell, qf_cell, uf_cell), + fe_face_values_ref(fe_cell, qf_face, uf_face), + Nx(qf_cell.size(), + std::vector(fe_cell.dofs_per_cell)), + symm_grad_Nx(qf_cell.size(), + std::vector > + (fe_cell.dofs_per_cell)) + {} + + ScratchData_RHS(const ScratchData_RHS &rhs) + : + fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()), + fe_face_values_ref(rhs.fe_face_values_ref.get_fe(), + rhs.fe_face_values_ref.get_quadrature(), + rhs.fe_face_values_ref.get_update_flags()), + Nx(rhs.Nx), + symm_grad_Nx(rhs.symm_grad_Nx) + {} + + void reset() + { + const unsigned int n_q_points = Nx.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + Assert( Nx[q_point].size() == n_dofs_per_cell, ExcInternalError()); + Assert( symm_grad_Nx[q_point].size() == n_dofs_per_cell, + ExcInternalError()); + for (unsigned int k = 0; k < n_dofs_per_cell; ++k) + { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + + }; + + +// And finally we define the structures to assist with updating the quadrature +// point information. We do not need the PerTaskData object (since there is +// nothing to store here) but must define one nonetheless. +// Note that this is because for the operation that we have here -- updating +// the data on quadrature points -- the operation is purely local: +// the things we do on every cell get consumed on every cell, without +// any global aggregation operation as is usually the case when using the +// WorkStream class. The fact that we still have to define a per-task data +// structure points to the fact that the WorkStream class may be ill-suited to +// this operation (we could, in principle simply create a new task using +// Threads::new_task for each cell) but there is not much harm done to doing +// it this way anyway. +// Furthermore, should there be different material models associated with a +// quadrature point, requiring varying levels of computational expense, then +// the method used here could be advantageous. + template + struct Solid::PerTaskData_UQPH + { + void reset() + {} + }; + + +// The ScratchData object will be used to store an alias for the solution +// vector so that we don't have to copy this large data structure. We then +// define a number of vectors to extract the solution values and gradients at +// the quadrature points. + template + struct Solid::ScratchData_UQPH + { + const BlockVector &solution_total; + + std::vector > solution_grads_u_total; + + FEValues fe_values_ref; + + ScratchData_UQPH(const FiniteElement &fe_cell, + const QGauss &qf_cell, + const UpdateFlags uf_cell, + const BlockVector &solution_total) + : + solution_total(solution_total), + solution_grads_u_total(qf_cell.size()), + fe_values_ref(fe_cell, qf_cell, uf_cell) + {} + + ScratchData_UQPH(const ScratchData_UQPH &rhs) + : + solution_total(rhs.solution_total), + solution_grads_u_total(rhs.solution_grads_u_total), + fe_values_ref(rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_values_ref.get_update_flags()) + {} + + void reset() + { + const unsigned int n_q_points = solution_grads_u_total.size(); + for (unsigned int q = 0; q < n_q_points; ++q) + { + solution_grads_u_total[q] = 0.0; + } + } + }; + + +// @sect4{Solid::make_grid} + +// On to the first of the private member functions. Here we create the +// triangulation of the domain, for which we choose a scaled an anisotripically +// discretised rectangle which is subsequently transformed into the correct +// of the Cook cantilever. Each relevant boundary face is then given a boundary +// ID number. +// +// We then determine the volume of the reference configuration and print it +// for comparison. + +template +Point grid_y_transform (const Point &pt_in) +{ + const double &x = pt_in[0]; + const double &y = pt_in[1]; + + const double y_upper = 44.0 + (16.0/48.0)*x; // Line defining upper edge of beam + const double y_lower = 0.0 + (44.0/48.0)*x; // Line defining lower edge of beam + const double theta = y/44.0; // Fraction of height along left side of beam + const double y_transform = (1-theta)*y_lower + theta*y_upper; // Final transformation + + Point pt_out = pt_in; + pt_out[1] = y_transform; + + return pt_out; +} + + template + void Solid::make_grid() + { + // Divide the beam, but only along the x- and y-coordinate directions + std::vector< unsigned int > repetitions(dim, parameters.elements_per_edge); + // Only allow one element through the thickness + // (modelling a plane strain condition) + if (dim == 3) + repetitions[dim-1] = 1; + + GridGenerator::subdivided_hyper_rectangle(triangulation, + repetitions, + Point(0.0, 0.0, -0.5), + Point(48.0, 44.0, 0.5)); + + // Since we wish to apply a Neumann BC to the right-hand surface, we + // must find the cell faces in this part of the domain and mark them with + // a distinct boundary ID number. The faces we are looking for are on the + // +x surface and will get boundary ID 11. + // Dirichlet boundaries exist on the left-hand face of the beam (this fixed + // boundary will get ID 1) and on the +Z and -Z faces (which correspond to + // ID 2 and we will use to impose the plane strain condition) + const double tol_boundary = 1e-6; + typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(), endc = triangulation.end(); + for (; cell != endc; ++cell) + for (unsigned int face = 0; + face < GeometryInfo::faces_per_cell; ++face) + if (cell->face(face)->at_boundary() == true) + { + if (std::abs(cell->face(face)->center()[0] - 0.0) < tol_boundary) + cell->face(face)->set_boundary_indicator(1); // -X faces + else if (std::abs(cell->face(face)->center()[0] - 48.0) < tol_boundary) + cell->face(face)->set_boundary_indicator(11); // +X faces + else if (std::abs(std::abs(cell->face(face)->center()[0]) - 0.5) < tol_boundary) + cell->face(face)->set_boundary_indicator(2); // +Z and -Z faces + } + + // Transform the hyper-rectangle into the beam shape + GridTools::transform(&grid_y_transform, triangulation); + + GridTools::scale(parameters.scale, triangulation); + + vol_reference = GridTools::volume(triangulation); + vol_current = vol_reference; + std::cout << "Grid:\n\t Reference volume: " << vol_reference << std::endl; + + } + + +// @sect4{Solid::system_setup} + +// Next we describe how the FE system is setup. We first determine the number +// of components per block. Since the displacement is a vector component, the +// first dim components belong to it. + template + void Solid::system_setup() + { + timer.enter_subsection("Setup system"); + + std::vector block_component(n_components, u_dof); // Displacement + + // The DOF handler is then initialised and we renumber the grid in an + // efficient manner. We also record the number of DOFs per block. + dof_handler_ref.distribute_dofs(fe); + DoFRenumbering::Cuthill_McKee(dof_handler_ref); + DoFRenumbering::component_wise(dof_handler_ref, block_component); + DoFTools::count_dofs_per_block(dof_handler_ref, dofs_per_block, + block_component); + + std::cout << "Triangulation:" + << "\n\t Number of active cells: " << triangulation.n_active_cells() + << "\n\t Number of degrees of freedom: " << dof_handler_ref.n_dofs() + << std::endl; + + // Setup the sparsity pattern and tangent matrix + tangent_matrix.clear(); + { + const types::global_dof_index n_dofs_u = dofs_per_block[u_dof]; + + BlockCompressedSimpleSparsityPattern csp(n_blocks, n_blocks); + + csp.block(u_dof, u_dof).reinit(n_dofs_u, n_dofs_u); + csp.collect_sizes(); + + // Naturally, for a one-field vector-valued problem, all of the + // components of the system are coupled. + Table<2, DoFTools::Coupling> coupling(n_components, n_components); + for (unsigned int ii = 0; ii < n_components; ++ii) + for (unsigned int jj = 0; jj < n_components; ++jj) + coupling[ii][jj] = DoFTools::always; + DoFTools::make_sparsity_pattern(dof_handler_ref, + coupling, + csp, + constraints, + false); + sparsity_pattern.copy_from(csp); + } + + tangent_matrix.reinit(sparsity_pattern); + + // We then set up storage vectors + system_rhs.reinit(dofs_per_block); + system_rhs.collect_sizes(); + + solution_n.reinit(dofs_per_block); + solution_n.collect_sizes(); + + // ...and finally set up the quadrature + // point history: + setup_qph(); + + timer.leave_subsection(); + } + + +// @sect4{Solid::setup_qph} +// The method used to store quadrature information is already described in +// step-18 and step-44. Here we implement a similar setup for a SMP machine. +// +// Firstly the actual QPH data objects are created. This must be done only +// once the grid is refined to its finest level. + template + void Solid::setup_qph() + { + std::cout << " Setting up quadrature point data..." << std::endl; + + { + triangulation.clear_user_data(); + { + std::vector > tmp; + tmp.swap(quadrature_point_history); + } + + quadrature_point_history + .resize(triangulation.n_active_cells() * n_q_points); + + unsigned int history_index = 0; + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); + ++cell) + { + cell->set_user_pointer(&quadrature_point_history[history_index]); + history_index += n_q_points; + } + + Assert(history_index == quadrature_point_history.size(), + ExcInternalError()); + } + + // Next we setup the initial quadrature + // point data: + for (typename Triangulation::active_cell_iterator cell = + triangulation.begin_active(); cell != triangulation.end(); ++cell) + { + PointHistory *lqph = + reinterpret_cast*>(cell->user_pointer()); + + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point].setup_lqp(parameters); + } + } + +// @sect4{Solid::update_qph_incremental} +// As the update of QP information occurs frequently and involves a number of +// expensive operations, we define a multithreaded approach to distributing +// the task across a number of CPU cores. +// +// To start this, we first we need to obtain the total solution as it stands +// at this Newton increment and then create the initial copy of the scratch and +// copy data objects: + template + void Solid::update_qph_incremental(const BlockVector &solution_delta) + { + timer.enter_subsection("Update QPH data"); + std::cout << " UQPH " << std::flush; + + const BlockVector solution_total(get_total_solution(solution_delta)); + + const UpdateFlags uf_UQPH(update_gradients); + PerTaskData_UQPH per_task_data_UQPH; + ScratchData_UQPH scratch_data_UQPH(fe, qf_cell, uf_UQPH, solution_total); + + // We then pass them and the one-cell update function to the WorkStream to + // be processed: + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::update_qph_incremental_one_cell, + &Solid::copy_local_to_global_UQPH, + scratch_data_UQPH, + per_task_data_UQPH); + + timer.leave_subsection(); + } + + +// Now we describe how we extract data from the solution vector and pass it +// along to each QP storage object for processing. + template + void + Solid::update_qph_incremental_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_UQPH &scratch, + PerTaskData_UQPH &/*data*/) + { + PointHistory *lqph = + reinterpret_cast*>(cell->user_pointer()); + + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph <= &quadrature_point_history.back(), ExcInternalError()); + + Assert(scratch.solution_grads_u_total.size() == n_q_points, + ExcInternalError()); + + scratch.reset(); + + // We first need to find the solution gradients at quadrature points + // inside the current cell and then we update each local QP using the + // displacement gradient: + scratch.fe_values_ref.reinit(cell); + scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total, + scratch.solution_grads_u_total); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point].update_values(scratch.solution_grads_u_total[q_point]); + } + + +// @sect4{Solid::solve_nonlinear_timestep} + +// The next function is the driver method for the Newton-Raphson scheme. At +// its top we create a new vector to store the current Newton update step, +// reset the error storage objects and print solver header. + template + void + Solid::solve_nonlinear_timestep(BlockVector &solution_delta) + { + std::cout << std::endl << "Timestep " << time.get_timestep() << " @ " + << time.current() << "s" << std::endl; + + BlockVector newton_update(dofs_per_block); + + error_residual.reset(); + error_residual_0.reset(); + error_residual_norm.reset(); + error_update.reset(); + error_update_0.reset(); + error_update_norm.reset(); + + print_conv_header(); + + // We now perform a number of Newton iterations to iteratively solve the + // nonlinear problem. Since the problem is fully nonlinear and we are + // using a full Newton method, the data stored in the tangent matrix and + // right-hand side vector is not reusable and must be cleared at each + // Newton step. We then initially build the right-hand side vector to + // check for convergence (and store this value in the first iteration). + // The unconstrained DOFs of the rhs vector hold the out-of-balance + // forces. The building is done before assembling the system matrix as the + // latter is an expensive operation and we can potentially avoid an extra + // assembly process by not assembling the tangent matrix when convergence + // is attained. + unsigned int newton_iteration = 0; + for (; newton_iteration < parameters.max_iterations_NR; + ++newton_iteration) + { + std::cout << " " << std::setw(2) << newton_iteration << " " << std::flush; + + tangent_matrix = 0.0; + system_rhs = 0.0; + + assemble_system_rhs(); + get_error_residual(error_residual); + + if (newton_iteration == 0) + error_residual_0 = error_residual; + + // We can now determine the normalised residual error and check for + // solution convergence: + error_residual_norm = error_residual; + error_residual_norm.normalise(error_residual_0); + + if (newton_iteration > 0 && error_update_norm.u <= parameters.tol_u + && error_residual_norm.u <= parameters.tol_f) + { + std::cout << " CONVERGED! " << std::endl; + print_conv_footer(); + + break; + } + + // If we have decided that we want to continue with the iteration, we + // assemble the tangent, make and impose the Dirichlet constraints, + // and do the solve of the linearized system: + assemble_system_tangent(); + make_constraints(newton_iteration); + constraints.condense(tangent_matrix, system_rhs); + + const std::pair + lin_solver_output = solve_linear_system(newton_update); + + get_error_update(newton_update, error_update); + if (newton_iteration == 0) + error_update_0 = error_update; + + // We can now determine the normalised Newton update error, and + // perform the actual update of the solution increment for the current + // time step, update all quadrature point information pertaining to + // this new displacement and stress state and continue iterating: + error_update_norm = error_update; + error_update_norm.normalise(error_update_0); + + solution_delta += newton_update; + update_qph_incremental(solution_delta); + + std::cout << " | " << std::fixed << std::setprecision(3) << std::setw(7) + << std::scientific << lin_solver_output.first << " " + << lin_solver_output.second << " " << error_residual_norm.norm + << " " << error_residual_norm.u << " " + << " " << error_update_norm.norm << " " << error_update_norm.u + << " " << std::endl; + } + + // At the end, if it turns out that we have in fact done more iterations + // than the parameter file allowed, we raise an exception that can be + // caught in the main() function. The call AssertThrow(condition, + // exc_object) is in essence equivalent to if (!cond) throw + // exc_object; but the former form fills certain fields in the + // exception object that identify the location (filename and line number) + // where the exception was raised to make it simpler to identify where the + // problem happened. + AssertThrow (newton_iteration <= parameters.max_iterations_NR, + ExcMessage("No convergence in nonlinear solver!")); + } + + +// @sect4{Solid::print_conv_header, Solid::print_conv_footer and Solid::print_vertical_tip_displacement} + +// This program prints out data in a nice table that is updated +// on a per-iteration basis. The next two functions set up the table +// header and footer: + template + void Solid::print_conv_header() + { + static const unsigned int l_width = 102; + + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + + std::cout << " SOLVER STEP " + << " | LIN_IT LIN_RES RES_NORM " + << " RES_U NU_NORM " + << " NU_U " << std::endl; + + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + } + + + + template + void Solid::print_conv_footer() + { + static const unsigned int l_width = 102; + + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + + std::cout << "Relative errors:" << std::endl + << "Displacement:\t" << error_update.u / error_update_0.u << std::endl + << "Force: \t\t" << error_residual.u / error_residual_0.u << std::endl + << "v / V_0:\t" << vol_current << " / " << vol_reference + << std::endl; + } + +// At the end we also output the result that can be compared to that found in +// the literature, namely the displacement at the upper right corner of the +// beam. + template + void Solid::print_vertical_tip_displacement() + { + static const unsigned int l_width = 102; + + for (unsigned int i = 0; i < l_width; ++i) + std::cout << "_"; + std::cout << std::endl; + + const Point soln_pt (48.0*parameters.scale,60.0*parameters.scale,0.5*parameters.scale); + double vertical_tip_displacement = 0.0; + double vertical_tip_displacement_check = 0.0; + + typename DoFHandler::active_cell_iterator cell = + dof_handler_ref.begin_active(), endc = dof_handler_ref.end(); + for (; cell != endc; ++cell) + { + // if (cell->point_inside(soln_pt) == true) + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (cell->vertex(v).distance(soln_pt) < 1e-6) + { + // Create some object to help us extract the solution value at + // the desired point + const Quadrature soln_qrule (soln_pt); + AssertThrow(soln_qrule.size() == 1, ExcInternalError()); + FEValues fe_values_soln (fe, soln_qrule, update_values); + fe_values_soln.reinit(cell); + + // Extract y-component of solution at given point + std::vector< Tensor<1,dim> > soln_values (soln_qrule.size()); + fe_values_soln[u_fe].get_function_values(solution_n, + soln_values); + vertical_tip_displacement = soln_values[0][u_dof+1]; + + // Sanity Check + for (unsigned int v=0; v::vertices_per_cell; ++v) + if (cell->vertex(v).distance(soln_pt) < 1e-6) + vertical_tip_displacement_check = solution_n(cell->vertex_dof_index(v,u_dof+1)); + + break; + } + } + AssertThrow(vertical_tip_displacement > 0.0, ExcMessage("Found no cell with point inside!")) + + std::cout << "Vertical tip displacement: " << vertical_tip_displacement + << "\t Check: " << vertical_tip_displacement_check + << std::endl; + } + + +// @sect4{Solid::get_error_residual} + +// Determine the true residual error for the problem. That is, determine the +// error in the residual for the unconstrained degrees of freedom. Note that to +// do so, we need to ignore constrained DOFs by setting the residual in these +// vector components to zero. + template + void Solid::get_error_residual(Errors &error_residual) + { + BlockVector error_res(dofs_per_block); + + for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i) + if (!constraints.is_constrained(i)) + error_res(i) = system_rhs(i); + + error_residual.norm = error_res.l2_norm(); + error_residual.u = error_res.block(u_dof).l2_norm(); + } + + +// @sect4{Solid::get_error_udpate} + +// Determine the true Newton update error for the problem + template + void Solid::get_error_update(const BlockVector &newton_update, + Errors &error_update) + { + BlockVector error_ud(dofs_per_block); + for (unsigned int i = 0; i < dof_handler_ref.n_dofs(); ++i) + if (!constraints.is_constrained(i)) + error_ud(i) = newton_update(i); + + error_update.norm = error_ud.l2_norm(); + error_update.u = error_ud.block(u_dof).l2_norm(); + } + + + +// @sect4{Solid::get_total_solution} + +// This function provides the total solution, which is valid at any Newton step. +// This is required as, to reduce computational error, the total solution is +// only updated at the end of the timestep. + template + BlockVector + Solid::get_total_solution(const BlockVector &solution_delta) const + { + BlockVector solution_total(solution_n); + solution_total += solution_delta; + return solution_total; + } + + +// @sect4{Solid::assemble_system_tangent} + +// Since we use TBB for assembly, we simply setup a copy of the +// data structures required for the process and pass them, along +// with the memory addresses of the assembly functions to the +// WorkStream object for processing. Note that we must ensure that +// the matrix is reset before any assembly operations can occur. + template + void Solid::assemble_system_tangent() + { + timer.enter_subsection("Assemble tangent matrix"); + std::cout << " ASM_K " << std::flush; + + tangent_matrix = 0.0; + + const UpdateFlags uf_cell(update_gradients | + update_JxW_values); + + PerTaskData_K per_task_data(dofs_per_cell); + ScratchData_K scratch_data(fe, qf_cell, uf_cell); + + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_system_tangent_one_cell, + &Solid::copy_local_to_global_K, + scratch_data, + per_task_data); + + timer.leave_subsection(); + } + +// This function adds the local contribution to the system matrix. +// Note that we choose not to use the constraint matrix to do the +// job for us because the tangent matrix and residual processes have +// been split up into two separate functions. + template + void Solid::copy_local_to_global_K(const PerTaskData_K &data) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = 0; j < dofs_per_cell; ++j) + tangent_matrix.add(data.local_dof_indices[i], + data.local_dof_indices[j], + data.cell_matrix(i, j)); + } + +// Of course, we still have to define how we assemble the tangent matrix +// contribution for a single cell. We first need to reset and initialise some +// of the scratch data structures and retrieve some basic information +// regarding the DOF numbering on this cell. We can precalculate the cell +// shape function gradients. Note that the shape function gradients +// are defined with regard to the current configuration. That is +// $\textrm{grad}\ \boldsymbol{\varphi} = \textrm{Grad}\ \boldsymbol{\varphi} +// \ \mathbf{F}^{-1}$. + template + void + Solid::assemble_system_tangent_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_K &scratch, + PerTaskData_K &data) + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + PointHistory *lqph = + reinterpret_cast*>(cell->user_pointer()); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + + if (k_group == u_dof) + { + scratch.grad_Nx[q_point][k] = scratch.fe_values_ref[u_fe].gradient(k, q_point) + * F_inv; + scratch.symm_grad_Nx[q_point][k] = symmetrize(scratch.grad_Nx[q_point][k]); + } + else + Assert(k_group <= u_dof, ExcInternalError()); + } + } + + // Now we build the local cell stiffness matrix. Since the global and + // local system matrices are symmetric, we can exploit this property by + // building only the lower half of the local matrix and copying the values + // to the upper half. + // + // In doing so, we first extract some configuration dependent variables + // from our QPH history objects for the current quadrature point. + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> tau = lqph[q_point].get_tau(); + const SymmetricTensor<4, dim> Jc = lqph[q_point].get_Jc(); + const double det_F = lqph[q_point].get_det_F(); + + // Next we define some aliases to make the assembly process easier to + // follow + const std::vector + &N = scratch.Nx[q_point]; + const std::vector > + &symm_grad_Nx = scratch.symm_grad_Nx[q_point]; + const std::vector > + &grad_Nx = scratch.grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int component_i = fe.system_to_component_index(i).first; + const unsigned int i_group = fe.system_to_base_index(i).first.first; + + for (unsigned int j = 0; j <= i; ++j) + { + const unsigned int component_j = fe.system_to_component_index(j).first; + const unsigned int j_group = fe.system_to_base_index(j).first.first; + + // This is the $\mathsf{\mathbf{k}}_{\mathbf{u} \mathbf{u}}$ + // contribution. It comprises a material contribution, and a + // geometrical stress contribution which is only added along + // the local matrix diagonals: + if ((i_group == j_group) && (i_group == u_dof)) + { + data.cell_matrix(i, j) += symm_grad_Nx[i] * Jc // The material contribution: + * symm_grad_Nx[j] * JxW; + if (component_i == component_j) // geometrical stress contribution + data.cell_matrix(i, j) += grad_Nx[i][component_i] * tau + * grad_Nx[j][component_j] * JxW; + } + else + Assert((i_group <= u_dof) && (j_group <= u_dof), + ExcInternalError()); + } + } + } + + // Finally, we need to copy the lower half of the local matrix into the + // upper half: + for (unsigned int i = 0; i < dofs_per_cell; ++i) + for (unsigned int j = i + 1; j < dofs_per_cell; ++j) + data.cell_matrix(i, j) = data.cell_matrix(j, i); + } + +// @sect4{Solid::assemble_system_rhs} +// The assembly of the right-hand side process is similar to the +// tangent matrix, so we will not describe it in too much detail. +// Note that since we are describing a problem with Neumann BCs, +// we will need the face normals and so must specify this in the +// update flags. + template + void Solid::assemble_system_rhs() + { + timer.enter_subsection("Assemble system right-hand side"); + std::cout << " ASM_R " << std::flush; + + system_rhs = 0.0; + + const UpdateFlags uf_cell(update_values | + update_gradients | + update_JxW_values); + const UpdateFlags uf_face(update_values | + update_JxW_values); + + PerTaskData_RHS per_task_data(dofs_per_cell); + ScratchData_RHS scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face); + + WorkStream::run(dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_system_rhs_one_cell, + &Solid::copy_local_to_global_rhs, + scratch_data, + per_task_data); + + timer.leave_subsection(); + } + + + + template + void Solid::copy_local_to_global_rhs(const PerTaskData_RHS &data) + { + for (unsigned int i = 0; i < dofs_per_cell; ++i) + system_rhs(data.local_dof_indices[i]) += data.cell_rhs(i); + } + + + + template + void + Solid::assemble_system_rhs_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_RHS &scratch, + PerTaskData_RHS &data) + { + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + PointHistory *lqph = + reinterpret_cast*>(cell->user_pointer()); + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> F_inv = lqph[q_point].get_F_inv(); + + for (unsigned int k = 0; k < dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + + if (k_group == u_dof) + scratch.symm_grad_Nx[q_point][k] + = symmetrize(scratch.fe_values_ref[u_fe].gradient(k, q_point) + * F_inv); + else + Assert(k_group <= u_dof, ExcInternalError()); + } + } + + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const SymmetricTensor<2, dim> tau = lqph[q_point].get_tau(); + + const std::vector + &N = scratch.Nx[q_point]; + const std::vector > + &symm_grad_Nx = scratch.symm_grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + + // We first compute the contributions + // from the internal forces. Note, by + // definition of the rhs as the negative + // of the residual, these contributions + // are subtracted. + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int i_group = fe.system_to_base_index(i).first.first; + + if (i_group == u_dof) + data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + else + Assert(i_group <= u_dof, ExcInternalError()); + } + } + + // Next we assemble the Neumann contribution. We first check to see it the + // cell face exists on a boundary on which a traction is applied and add + // the contribution if this is the case. + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() == true + && cell->face(face)->boundary_indicator() == 11) + { + scratch.fe_face_values_ref.reinit(cell, face); + + for (unsigned int f_q_point = 0; f_q_point < n_q_points_f; + ++f_q_point) + { + // We specify the traction in reference configuration. + // For this problem, a defined total vertical force is applied + // in the reference configuration. + // The direction of the applied traction is assumed not to + // evolve with the deformation of the domain. + + // Note that the contributions to the right hand side vector we + // compute here only exist in the displacement components of the + // vector. + const double time_ramp = (time.current() / time.end()); + const double magnitude = (1.0/(16.0*parameters.scale*1.0*parameters.scale))*time_ramp; // (Total force) / (RHS surface area) + static const Tensor<1, dim> dir ({0.0,1.0,0.0}); + const Tensor<1, dim> traction = magnitude*dir; + + for (unsigned int i = 0; i < dofs_per_cell; ++i) + { + const unsigned int i_group = + fe.system_to_base_index(i).first.first; + + if (i_group == u_dof) + { + const unsigned int component_i = + fe.system_to_component_index(i).first; + const double Ni = + scratch.fe_face_values_ref.shape_value(i, + f_q_point); + const double JxW = scratch.fe_face_values_ref.JxW( + f_q_point); + + data.cell_rhs(i) += (Ni * traction[component_i]) + * JxW; + } + } + } + } + } + +// @sect4{Solid::make_constraints} +// The constraints for this problem are simple to describe. +// However, since we are dealing with an iterative Newton method, +// it should be noted that any displacement constraints should only +// be specified at the zeroth iteration and subsequently no +// additional contributions are to be made since the constraints +// are already exactly satisfied. + template + void Solid::make_constraints(const int &it_nr) + { + std::cout << " CST " << std::flush; + + // Since the constraints are different at different Newton iterations, we + // need to clear the constraints matrix and completely rebuild + // it. However, after the first iteration, the constraints remain the same + // and we can simply skip the rebuilding step if we do not clear it. + if (it_nr > 1) + return; + constraints.clear(); + const bool apply_dirichlet_bc = (it_nr == 0); + + // The boundary conditions for the indentation problem are as follows: On + // the -x, -y and -z faces (ID's 0,2,4) we set up a symmetry condition to + // allow only planar movement while the +x and +y faces (ID's 1,3) are + // traction free. In this contrived problem, part of the +z face (ID 5) is + // set to have no motion in the x- and y-component. Finally, as described + // earlier, the other part of the +z face has an the applied pressure but + // is also constrained in the x- and y-directions. + // + // In the following, we will have to tell the function interpolation + // boundary values which components of the solution vector should be + // constrained (i.e., whether it's the x-, y-, z-displacements or + // combinations thereof). This is done using ComponentMask objects (see + // @ref GlossComponentMask) which we can get from the finite element if we + // provide it with an extractor object for the component we wish to + // select. To this end we first set up such extractor objects and later + // use it when generating the relevant component masks: + const FEValuesExtractors::Scalar x_displacement(0); + const FEValuesExtractors::Scalar y_displacement(1); + const FEValuesExtractors::Scalar z_displacement(2); + + // Fixed left hand side of the beam + { + const int boundary_id = 1; + + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement) | + fe.component_mask(y_displacement) | + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement) | + fe.component_mask(y_displacement) | + fe.component_mask(z_displacement)); + } + + // Zero Z-displacement through thickness direction + // This corresponds to a plane strain condition being imposed on the beam + { + const int boundary_id = 2; + + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler_ref, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + + constraints.close(); + } + +// @sect4{Solid::solve_linear_system} +// As the system is composed of a single block, defining a solution scheme +// for the linear problem is straight-forward. + template + std::pair + Solid::solve_linear_system(BlockVector &newton_update) + { + BlockVector A(dofs_per_block); + BlockVector B(dofs_per_block); + + unsigned int lin_it = 0; + double lin_res = 0.0; + + // We solve for the incremental displacement $d\mathbf{u}$. + { + timer.enter_subsection("Linear solver"); + std::cout << " SLV " << std::flush; + if (parameters.type_lin == "CG") + { + const int solver_its = tangent_matrix.block(u_dof, u_dof).m() + * parameters.max_iterations_lin; + const double tol_sol = parameters.tol_lin + * system_rhs.block(u_dof).l2_norm(); + + SolverControl solver_control(solver_its, tol_sol); + + GrowingVectorMemory > GVM; + SolverCG > solver_CG(solver_control, GVM); + + // We've chosen by default a SSOR preconditioner as it appears to + // provide the fastest solver convergence characteristics for this + // problem on a single-thread machine. However, for multicore + // computing, the Jacobi preconditioner which is multithreaded may + // converge quicker for larger linear systems. + PreconditionSelector, Vector > + preconditioner (parameters.preconditioner_type, + parameters.preconditioner_relaxation); + preconditioner.use_matrix(tangent_matrix.block(u_dof, u_dof)); + + solver_CG.solve(tangent_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + system_rhs.block(u_dof), + preconditioner); + + lin_it = solver_control.last_step(); + lin_res = solver_control.last_value(); + } + else if (parameters.type_lin == "Direct") + { + // Otherwise if the problem is small + // enough, a direct solver can be + // utilised. + SparseDirectUMFPACK A_direct; + A_direct.initialize(tangent_matrix.block(u_dof, u_dof)); + A_direct.vmult(newton_update.block(u_dof), system_rhs.block(u_dof)); + + lin_it = 1; + lin_res = 0.0; + } + else + Assert (false, ExcMessage("Linear solver type not implemented")); + + timer.leave_subsection(); + } + + // Now that we have the displacement update, distribute the constraints + // back to the Newton update: + constraints.distribute(newton_update); + + return std::make_pair(lin_it, lin_res); + } + +// @sect4{Solid::output_results} +// Here we present how the results are written to file to be viewed +// using ParaView or Visit. The method is similar to that shown in the +// tutorials so will not be discussed in detail. + template + void Solid::output_results() const + { + DataOut data_out; + std::vector + data_component_interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + + std::vector solution_name(dim, "displacement"); + + data_out.attach_dof_handler(dof_handler_ref); + data_out.add_data_vector(solution_n, + solution_name, + DataOut::type_dof_data, + data_component_interpretation); + + // Since we are dealing with a large deformation problem, it would be nice + // to display the result on a displaced grid! The MappingQEulerian class + // linked with the DataOut class provides an interface through which this + // can be achieved without physically moving the grid points in the + // Triangulation object ourselves. We first need to copy the solution to + // a temporary vector and then create the Eulerian mapping. We also + // specify the polynomial degree to the DataOut object in order to produce + // a more refined output data set when higher order polynomials are used. + Vector soln(solution_n.size()); + for (unsigned int i = 0; i < soln.size(); ++i) + soln(i) = solution_n(i); + MappingQEulerian q_mapping(degree, soln, dof_handler_ref); + data_out.build_patches(q_mapping, degree); + + std::ostringstream filename; + filename << "solution-" << time.get_timestep() << ".vtk"; + + std::ofstream output(filename.str().c_str()); + data_out.write_vtk(output); + } + +} + + +// @sect3{Main function} +// Lastly we provide the main driver function which appears +// no different to the other tutorials. +int main (int argc, char *argv[]) +{ + using namespace dealii; + using namespace Cook_Membrane; + + Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv, dealii::numbers::invalid_unsigned_int); + + try + { + deallog.depth_console(0); + + Solid<3> solid_3d("parameters.prm"); + solid_3d.run(); + } + catch (std::exception &exc) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Exception on processing: " << std::endl << exc.what() + << std::endl << "Aborting!" << std::endl + << "----------------------------------------------------" + << std::endl; + + return 1; + } + catch (...) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl << "Aborting!" + << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + + return 0; +} diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/author b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/author new file mode 100644 index 0000000..175761f --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/author @@ -0,0 +1,2 @@ +Jean-Paul Pelteret +Andrew McBride diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on new file mode 100644 index 0000000..f55642e --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/builds-on @@ -0,0 +1 @@ +step-18 step-44 \ No newline at end of file diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/dependencies b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/dependencies new file mode 100644 index 0000000..e69de29 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/disp_soln.png b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/disp_soln.png new file mode 100644 index 0000000000000000000000000000000000000000..3859e040e4ecc4cc2a0c6166538bea9d5efbfd8e GIT binary patch literal 28134 zcmXt91yEGq+g`duIs_M_JC?4cVX37-8kQ~zLAsHUZV+5bX+cRzk?!u24naXm%J1^~ z=0A)$v)n!BjpwcBoEr<*R3^Zs!3BXp1ga_uIv@}#8VH1DiH!<;B3*L_1^&adR#R30 z{d@dDwtq|mKEZjZV(bM1;gLT6LjmQ!qXIs}@>bPQ#9GHE#}TD*_{et+0)auQ3bJ~y zmyYxO+Mep?bZV`;r6wojqvK;oGC_^rQnh%>3NJD-x;1}O(CcEujwwE)rF=<63||Z< zNq2pV&lry#j`kg+kAaD}7n_Q~iYPgOI1SoZW6~+Z;dtXS?>YD4Yr)oqnA3d2t<|IL6z2f+;pt`UlQQDt$)$_$un zD4_JGJs@UKr~pm)xGNue(bvMsL4Hz%{r?}!rM_W)69w_JNeW4$g0T6~?F|(4vE7S7 z*bE%BA!OCopx3A%w)d=E6o5bzKcY8GGJ=0dbKNwRbypFA{}bvP=}b41)s|HZ3Je6w zx5n;R{#qvLEhX#l!C_9}KY_oH z`u9%@0O8~@dCHgq!r>Cr|RAlyX?f}-GmBhQ#)uoiMWb9l?@0O9;^O*oT0Jb`}T zjkHO8&}^Ume{05BV%O9}qTj!$V(JM9&;0W%mTH^s5FNCW8FnPh15O^iP@BVHBiR6K zo2-}SsqnP`Wb$a)M1pVUu-7X|FN}!&)!$k|W=V)6B;|duK_V09?z{(~m*4j#!&pI~ z;R1g;If6%Cn)7NfbsOYZ;lWWs1I73&<~Y{%&z_w-lPtm0NQoFCiNE+L<3Ik5o5l^= zN2+I1ipoF0FoOdcvd2DgTdFq4wupf_^tgYEpTo@g4$=aF{8GXKIDaYp&h4R+v8%Pn zkgocd?KNlQ^5_b`B9doATjfgbFNgzP^neYM=}_6y7+xb^p~xyS?Rd#;wSM%h1=~gf z0=^sQlFt}GvC?POU@FmpR`<`>VgPvCPKs19)mME`U{~P!qVycICG;mTR1E*y+uAm` z1WiuSHRV(S;iD8ADb!Q)zbej}JJ#`>hj@D)otKPx(WpFA+P5Ka%+ zYk&~NTeMw}UCT$$$zVR^6dMkRt=^L^yYb z!)Wqqi_N4Ktr3V(5O>Y{6h)A}u`EjnzI8GS;tfbV(n4u3Dwpni^uB9z&s%gtz(ryT zETSgEJ0${JUUP&Q_vWX6hd+Q88A3HT>;(r(K?NXWk12mK<3EkYOz>1A6UV(q9R?u0 z%I?&I8`AG{@)$m&E5fN;Te(kDKDR;22Mcfpht27j<0AAA4iN2#W)2@w4G#6yQ^@B158?>d1h0 zLmC$MZB(iMEOA=>pl^N5UgCkn*gn~Gac$a5U4lh!6B<5g`ZQBa)H~LA$-$d7;WNf_ zHacKG*hmh^PT1EPHNJp^IIO|%7mrqmc}Kr$qv>9|?m`*JxlDFMjd79B`)G<#fwS30 zMZ{UJ6>n5w<64X4y>E_C(dTtWSU(#Il>Ic-+P##!k3-op?SWjAJj%095WwGdzOKC< zB2{K|k$!NWt^n+*Y`tpqz@+b~04zB5JBm7^%zzNqYc)&$t)-Swov9bdF0oJat#dYO zC+oGV}I4S z{0_}&IOO~jx6Pwq`EMz2%`~bQ@Y?n58b7!%dE%tUp8`4UM>dPR^XN#jI-L!+G4Uc^)>ucpEmIy&7?5?Hja`% zF6r5Zv*a?eMm%Vs7TD=B%VJGP z6$`D+=|}&KeSo@Sru^jl>{LzyPYbT@p|yHcRMFa6H%W}w|81pItcqv9P5FeYZ#J)Y z){`0L#40cQA2=Y~?nvs4T?W5>Q7K5Evn09rT@oeE=xBn_5QzWNtPs{0;r_>r)nfl^ zU#wXBZ%z+~N1&0GY3zrUO3hF4rn8oRjW)ffKG+bWF$3^HK=s`E<2Aqcl-SpMDd6HT z9{+^To>o>IvSOu!z+OJi%i&6Fb2#d8=mJOhYhtn6YeGxI(fR?2L)cS&Hg>M&i|%43 zl)!#}ZPa@-5wo?E{xbv36m&$+bv7YZj-7&|-fIeP$zx<)hE?FOxF-F zlRa_yhwD{3UR^PyODM1sw}1h(Toc(o(3N(CU~TnHIP69Zh#%+%#K!nu`xuMri*`z- zG{TfBpEWmn6C*OzL#OZ=8^|$z@>J%x2{k3FhF0GI3!#{dVSj4$cjbTKSXWnp!?%8@ zWA!VyVljFBWk)i$3{YBuhoVQjFW;XsjfBCCU8>BOY6W@+T4a*Nyf$(8AF-%TMN{|h z?a1>~qCS@w##)Cuth#b4g!?r5DDn4En_2Aj?av=kmCTVhRu)NU@2t8DE#nfN%T_3n zx-(p8?VB_Ztk?zFKkk%_;tCX7+-xpk;l&}PzGqCICQ>a_xH@??G@Qsix)*OOsQz1u z#N34%cy2(85w?FMUvAkRCvSb!o^qGhXOs1B)ALs$CiA6kL)n`w@xVc`{)u`S?e4W(}8A*wCAA;EcK0OrR zE3-YPyC|*&0Bs|qHh2?eFQC~ZaX8P1^<@0NC{09Is%)(X6-_bUWy`kPf70Ag3``(p zrJS9)|Asl7+U%>L;l||V%uqM7(8=+`5ZKm|1?_vW4(+i4vc{G$-Hw9&72Da{Jxb$C ztUj~oRH7~lpu2Gy?xMK=o2uCPs@wWQ& z-d{wXy_DFdV(?HG@)^h#WiMlykH&t)NG&Vc@|&+Hr8;4K^^ApU1u81PVp+C8YdNwA zTY}o;k}Skj>wLbJRhZ2*mVvS3RDsrl!|e< z0(wqBv9}0IRvZ!dF{yd69?}UzeMB}e1?q^Y9eCg6?xs<(`cL_l?q<(q`=PbXJCp(t z)C^Netr=xnRYAa+F#W0g%yv<>shbVlw3Goa=rMreT_;f68i}?~$$D$;>r<|ljmD57;UONu)*-nu zJpMQK+5Eiso6G23%dM@)xD9=a|C9X($VHQ}f|g4-@eD<%w?Arbn-TA!>*Nx_KRy)= z6cI2t+HKialv$myXt88tm!Is7tW)o$5etlZ1EYUqE!?>*-eU(;Y(HY0m>2C=!#-wP zEkn7LJ0pP?5^(CdKJAON!%@1=pKr4VO=eOx6M=ql=xwYgA$W@5TyeZU9SBfY!t^E5 z@-Z-MVZ~_eS@wGdyLOMFR}Q}vAhm6e3xnl_eBd;;WNWc&bF5DW-3^i3Xb3xEeKY)j z8H8oDrl9x)pfqtrBTcMlc=d&A%~|<=X_VBR8qc&nUzaV0JMV%|7isyO~cjTmy)9{@q7CIMJR#JV5-DinC8eS2xu6-2yZDg4mo! z%B5aY9#o6)R6k;}XE-z$ZNkL;srb~qi3^uIVPrzWILV7Nr@fOX565rY#y}{x*lr~6 zYgbmF!S~f@WHjfO$3_A%Q%25EZ_(~5&Xw$ZwAQ>And^bB5h{-<&Gh7ynRfhOoH>Ce ziBmC<6Q_s@C1yIb?Bi2zT++*KAVOzzY7#J?HS$R*~Q-$Hj~hQ z#J@Om(K=_A zhf7jn6fad0{DZGGEk<&cz7!@Hv)M|dc=$1+Mcyti{qLGYBa-p7BbeP4K=}|C+Ixri zr!cW+`G+1OKzSDr3k+vFjp9$D)p({8B+Pk4eA=<|*gU5jU?yfIo;6nVdw5v9dXl5A zimds(nlOs4tq!hPG$1_qTIaw^H{8Yx15LCD{5Fo|?+A?XmB@m4)YVM|GUQ~_f1&sm zFHBm{OU&ECA~yK+I&QbiJzza`0W}>2<3Wq{V~`NA6v%=|E6p=2-6xVLhemzG_BBG5{KQjT+we}=p!@5 z)YN(VP=OXAG;QfYqHWgXX~`)|Mx7XCoo;AIXo0<+RmnMXJ}TS4{3tOannfn|il!RI z=d>7n6ZmV#oQO#DlnS)kRiS|-scPQ0B7?Fxpg0z^kBw1VWH5}&z0|{3+_OhktU4;h z2SU+!c{JO%o=@5+mJ=B-gX|5&7Xf5UoC^HbZ^?FR(%V4rIJT)EQo&O>(j`r`U=Js_ zhFCR|#N8)cZSvKQBtTv>Oh{4ZrT-Q_%ag&#pFIOLDPR(-DG~rk@n=(J5xHouJUd$| z7FnDD!+9V!&P$N)`}hc4e0_cDGS<>0NN)BabIjeoY8kGX8bNX3xo(Hj^@)YZrM#|u zsi{#Y3H0MoSB|QzMTZ_z{I`~lHl%j{NmTW74;Z;iw6G+MI?uOIp*o?!1c#r9Dn^j!neC7S{S?6aF*#(C3*3A8qi&T z9DwR?KJy>znv)hzDQu&sGq?RZIoRX-VAk_Rnod<4 zb)Sn`RQc3o*>)HI+!-g;OO>E?GKVP!0ReXim1B4j^T%n8{nO(<_*B^2&X7=j8yfwc+2f2?q}M1d+Pdb~fbP^$EdN_XX662XiA{9Bpf3jKYs?7)mI} zvosi-(8<`4n5I+b85bo4VU(vy*WWkCDK37Mq1{Zfa+uM|zn>@bX?uwKagn6LB7+&$ zI!Dn^Ns<&(OBomQ+=TwTlD3l3Je$RTQef#keLWDvyz1J&?(0uxoU1;xuNJ3jXZhWv z&>x76k0iR@=}wq5rKk0Mg0C%__@`$K!UuF;c)=6$B`YASD*F`V3Vv4=3OLAp zftM}*Zp!i;I{6AEO&%J7Z8ghq3_pV8UshSk<0xlEuXZf2OVxDpN$D!?ZseB%6x zWqJ-1FBROmu&`Mq*eVzk`A+P&d=t1vNe_F5vUt&0f5IDjQNNHMhZNqhxt=Hf(OmD7L0S+8%h;b#u{ZO zug|AS-0NF_y~wQOZWWd)j7?jasEr!Rh1j4WN2%Sxhxlu=x!SV`WzF=kv0YRWf5l-e z#lYgvXGP@UQJ;j3bOTu?jvZ}2x_7jP@7UTetQVmzQv}0E%$;6@5wkmndhX%O_G{2g zMbS-2&GXJu5-^8K4tEiFew!I$D&qOF{|WM2Bv-8hK!M}Nr;ns^{Ic6k8HvxaZfCc9 z9fU{2EI-nOnhuBFv9%_r#W2~l>v(Ng>l7LdD&ytaFWAaDrtQCDAIQU8 znK{>Rb<`l9Wd=WZ3`&#A-B2-1>dQTZ2LHRP&sDF7B*Twoend6AhzTj7{}6CY$1W$K z6?qjI<3v_?)691y#;~u5o7|3nrxYamDdk_Rc!pAW*9XC5E#g7LBmz=}?!zV2(lOFS zE-tVHdU6B)3P%cF>IhkH >fV=AsBz}>XQkC6oA;%4}N?%w>I3ybQy><&Nazb6U4 zg`*B+;8Wp38N{DNJr~TOc1&;^|PL!i7By47&L+z}9DhYm0-TZwdGs~JJ zMMp%O;Y!Y`+%TVSc?9ZO`~Kw7PhdLgX5D$*Y%RuAQ+?zQ{Pa00)?uiDS?sw1*?j3< zMs&}*IwghozxNhnpBP}E#i6cMl0wrn<<>lkY(&D(apu|Nuy;%`j06hlm@49!KkTez zjjbry5v3?s!cj3i3^z=2U#3c?dX7iVKjh#Zn2qZZ!3)$58%BRtsY;UHRH@{?=9qc* zWGdjjtLQ*HwkE}1#@whRQvs-}5x?vO${<4jfen%rSRH^b-ZeeLeht;^GZec;qR}gI&htbcfR!7|K6D%Ec z#(+jN(Ca*3ips1G!x*Z>BAGM^rRWsxnsm!js#dLxGkKBA%LUMUQ^{d67fuQ54RxPS zxbhMbZFf_5S!KzDtn^f#elQhJ1YjI>-S`uJ_3U0CL_`}^#R9Wps>!-v6irAqIY%w` zhg%lS6ANKky@tB4MgfnNV#5HF{!ss!=+Qw~ET2l|BOc{?7jCR=81B2hNfgW-!c}Z) z*X2gJB=p>}6>$@VHp12xVsb^tb!@iXL|gT!+DsqjC*P~l2r|z&=cz58*>K-1_Usm) zYA6-cr?_eys>-nW)yCu-=3a9F*2~3TnO~;CY>T7iM2g3jRa_7a3L2?UP_sl$sT2Es zZDpL#Z)A%5JKV$^^EKpG)?H;B36w(sOpcu9B6pxA%=wee5|2+w`|%x^z&HeFKJha;%S2Rk1z{7?EWRs>Ql#-Cp1bU)t*eI#JdpzR~(W3f2eNV?c&M zfYE%fI!Irnrq$AvGbE7WZnem>jCFbYg=GjUbZ;onqNhBp3GJG^7@B#L7pQbqExm5c z(PycTcJhyFPBcY zc|Tn@jk)0~0KFRyow;r-lE5^)nba}|`Hri&CMJ)Mszg8PQ$x|UHM7LDK3|^*9 zZ6t1tLkE;T2+~0BQH+UE7MArpPe$%<2;S;3?j|KNwxpU6K_9y^h@{E9@@;Cy*Y6??6{3ds ze%U0_Q$DF+FHX;w6>30F8ITQE$!3aS8Vxhui#<6~eH4ieZCZo^^AL7`+^B_l$J94> zsGCqBU>Gz)Xmt6t3khRh!k-`!=qbG#cmVUKLET&@;wFO)WkRsR@dsd1>}($#5X;Ziu{IZ^Oq=)~g(RPMYyciM7SL>m`4f z&mvNYs+a3|*&C^S#h%9fGw9X03c(K?s0rqItxU#1&r9rpQq)1x!&l{%N(rs{X>zH? z`RXK(Jdg4&2#yL*4zoC464N%-oyuMo7%xR_y3O@FEfl<@Y6nzhlx14Zow44-f#N`w zF(scoto!$D-PrhOm+6F&+Y<*{vBW1_^xanEA0X_o$}8kvBMgO5CcP3=2x5i(Fk)sy zD}^aAA>ALeQ@X?*GjW;L^{ht2M*4)wu zk56Jx+;HbmEyNx_{8Kg`v;^6hijmsSTgGbmC!zvslHnm2m$C=?s$g$D4n9X!vI~jm zr^v}@whtrVk~~VB40?rS~iW#|!x^@?M5?qIuallUwdoqYW8W=ZY1%7X>Mb?$w zd&6GJ4coSSh80Q@T8>|)__DPE9-|ZHG$C?x z4Akn_I`jQZ;(Cec!4bx4)&$55iNYC4IA1bdxF|t$S=0L#%^!GQ7#Q@P=*#$LPqAk4 zratTCbaiL<;sAEKtJjWdDWOqi_#0cINvfbtswD*c)>RboDcF~OKLhZ@MTm`Xb39}4 z^C@EsyAymw&2*ulim6N0WWt!>z8m*?zuW+Mn#O!qlynRcfB^$cM5}85Llgb2i}^Y9 z72W4v4cJF~C$??co50l!4~=ck`&$E=gM^i9$A(b!!jy&wVR}6D}D# z;wonXB=q$S3!>7pOxxlG71EVC;J6pURiC?>Z*NPYf1^(@qRAv1@oW@1a0S@;g=|1{ zw`Zsa^5}2l6PQ^kae56tT#_p2i#P+J6iOhl!K`@b4L7&2;~_drs*DaOH5_D(=VA&O z5fM-mlP#Ic)!S78B%+0Ew6^1TD6!))ed$o=Hk-D-=3HV=UwAgk`lNIsnP>Z3U!PEv zCNZ;c#g~c9n&Z=xf67h(MIvh?rEEQ{dc^+xe#%3UIix5okf5dFBB<^JzhZX`OoUT( z9oD&XAfFAKpvob@9{Q;n`|f)#XVq_Xc*yni9}1b|bcmdIg)`N$J^*^kfCXfkm`D`V z8JH!|?;#r|s?>!&cS#bEk0dP%>qgT5iVr30o1_W-YW?iDCWy3HWJ|JmJ@dhF@8Koy zz)t?Pw|S2yV$KTt@=8#Snk<1}o^-!7T5r63rcj*OeuMoh%HT@MrQN_HQ%uEFP#j1L z&PX8oMWiMmGRu0>=%Y}1qfG5o?o1}bgbyZ;4}%@q``4jA=O22tdfZMU?#boEv4Y1D~>f0pQ{Eg zn`{eW!}~(RPKSpE@6Q+g*(81oean~(4VP5IoLK?{lINbJ^x2ktn@hOLb6a%8o1rok zi(~!Yhp?PyzDy#_8RL}fJ1^^<2ND`d%U!LDTZ}bkPE&Ch3kxK4?k8Qt%Q(Jv$^kBc z%TcQYSi0p`Y^p3u3L)wQ^Ghu1+S*I^uGh<%V;~|T9|qA2K*A1P<*v6HL(M3Eis}wE zGbv#GOAlSEB~#fct*g=AKfF{ml5h_I;(IMAq#|z>`VcklZ&1hwR|Wm zTLgadz(NV6OGI@ZJL00R#Z^B)L6A?2xO6eKtY9A>T^evAJ%MPQ;pMN);K}A71*eY= zT%8R;xoNNr{4h>?Wc$Lx%U5+iL`9eP1A6leZ$e1O;!e>^(@;2Bx=I{eoPWAoWo?;e z;HC%T;f+B9>gq~tCXTET%{^R;oGLDED?zB`43e_SPF;I?H?dA1FODXCBagymyKe;6 z7i)wDw@!Y(3rlY6?n!AK91~{L;*KQ|7U=^>VsT?k>3oRRjIQ6DuB1HXdS%5N3Po0^ zql@)cH#OTb#?csb9D-db7;{cbj8cI7E`%(Ym&wo1=gL$5rfHG}BfLANK^-3JB4tP$ zko(8|N11QEuj>Z?ar|{PK;G&;21z~>+WVjS^hunf3$;6#+tc^vJt3vtMeZ9!p86cM z2lzKHuCDQ^TUZ}li%;;r4=OU2J4aQJC!3joUWXG-h9+IG)E~eCO@r@v>NTC4T)+Xs z`XNYr%6dcceW8*s(;;#Jn{d5y=&q!k|8w*d$f#;V9I-yrenBKd7peDK{DGt!0wTOo ziJYYS1ieB=sdS3V#2G$x_VcRDpt%d{t+Hi~K8V9O_>5_ePc4H0@RS~bXrrM_ZAtBv zsJ>i0ddIn?@ps&w0KEBVs2{WJ*0&i7T9aSKtH9XQgqo5^D`a(XKJqv6lz# zXHBAg*9jZ3%M@5bz7WJ-Y*uO|NQu_N7NdqcIdkY;nG~r&eaaRDUcPp({B_|yEt{J2 zZnGNVMt`=`B2<1y>!>ACQ6ZhTH3yK0}Zw0Prut%qBXVHg>E^3DQ{|oN|y+ zGbv|iS5FzCvS4>+60%TnR}IrK=$8Zv&6%PAPEk(-#wu1CO*bi-3MQe3?X+JjX9Dv& zPy1YnnH|83pV|DSP+gLWFCoM`R-+M@9hk#sl$r)tje*b5qpT%msUc9$*98wSo4sn) z;@nj6GprLBMujo*PCVk;qrJqDm?O$|@?Sk;5$a?R_%VNCQ;U|(pHT~1IKCtCS4170 z{CcNWuZE%jB}Aow!>ZY})Q&i5(?A*eMG$6K7aH@=Y0A|WIM7N7{+X$0zUxVsx4$&cZ4Wv=XgH3w% zm+GB}U;3);F3~$54(?PzV$9)*NYflyqt*N7c6sgbhsXH%I}d7I6>muL?Lh64b=J|y zLe?#*4e~VekPnDPG+V3;jrA2`sU;ri=E*9QQ)wmJ5bU2+dtZqMqK&WBu1*j%yJ;Ni zZ%42F{V8KLOuVnfiA(42tWKOty*vC)IRH~zo%{zF9|e}IMP2*8o_}P4tixgD=|NQ$ zukO+!lgAZo7TawtdCXxew+H%^gFN*h|D04Ilw?~SjpfpSrNRV7Yx9{8Ipcjs->HUb zf_0IERINV8%gJQTL$jZnPj%?#KMF@+!@~myLtl`MxlS%*HUS4>!#rOmK5c`Th|<dk9^aAf05Pm+ zkGyhk2Wad+vBiAoY_jwS@L!}gd4g201;f@{NzPGN--RCc{iY6M%96EBX$W2~dP)oh z&{GS{hDrsLBB0V|cIJ?|KPWWw&%dOVE{)knSMMXt9oqGhcTA^V$2O6ljFKwE2JI() zBSxYx4uyWW7k&zp;HC*6OeD+sWnb=!QMDiWN*=k34=`W8V_lBnDm@*<_Ms#QCyYx- zEh?P_q)|``az4UxG+cktrCoVwS|yMc!uJM;@g@2yODe?j%DnuO)lmGp*-ytRUqnS` z^IC9b>L>EXTP-7Zi6>sW%>)S(Pdq*U^!6yU#54F{Z4(c4R7SkT#BwoO*1*^BAzg}LP56p7!?(sSEFQZ&c9D{4 zGO9mkp^a!NR+vwa=wk>=?@te3xhYXCgA3H-`ZUXEp6)uNghdP1a^ib@{M@`C%8FAo zZu2Mv1Szm2E8^%|8s1h7(Zu<4(MO$$t7C4vJ4cPk0HdaNC1C!XFXli9K^`1*BZb6De5EJTt?RBZ3lJyWdC2a(qiyNteMIQhGVIj-gHBrtD zwx_y~F(vYfIIH#;fbxv<*c<~EGF26fEJX(lBFWZ<@T< zZQZ#azh;#jPHwBd1V_ z8Aq)(cH|KM=diUlJFwv4{8}$?9ZiscYSM&y(wTcRzDf;4jGe~h?E|iTnhPv9m)_wh zvd@$y&r4MZlkZip;UMC*3UMKNxpQbSic<7erb~QZ@7lPI^wVW-a)-zO2`4MSPm|Ot znkcE|98J$IU+KuRg>~S&Fdw5$g%#C#>A$S&P~W8j?!A39p;RH*YhcC8835rhrgvpv zG7W?l^EZFIs0d88nq1IK>2h<{Aw}8OkzC8GyqB-<(Uj&J3q z+&jOn1|7&;b_eU3Aw=`FtKyE3k4YK_#(AKd6ASZsxvL>K2S5b$1)A`$YL~UVsR|x1 zG}+08vD*1Ip>ov>BP9@Iv!ifKv0!n|+J^GhScy5NBbEV*la2rnCmSM*>iUMg%pe9j zQ2;xmNf0(6ZsgPUe#r-e3=@uOLqIi`2S(d{^&4=KsG&&KMSP*d3e7E$*wd0O0ESn`30}ef+xI^?R&f_-e!plj`QsyzRpan( zpIpF2KR=ftM$&8FR7c;Hn6Q!4;o+(U7|_FmY4~cJ1o--cO6XkAGGfJ&W$fwW`Y7hI z9P$4MF?JljUauj7FE{tq($YP8y$UQrdQn}*g{?sUwx3RUhhXyh?eQSK# z!3?C%QT2@f7iub_@v)C>9U~`77zV+`=X9w%NcZ-Ck@1AHN=^)7Su@jL@qmN?*I|CG zjE4U!{yN4WJ*7FmsO|#j=@u6SIB6s}$Q^hE8?CV+l(e+F;q{~nGT(^Pb#XW{gbCod zUQyYOX?sO=c$~%Ow|2S^pRe&bYrfp~s--$1JnaRf%_!+mb2R!gJ^BGgL*(_&%b$rJn7*thx?HwDOOW#zYbWx zfHDV8=%rm|rY{%V%;hIU_FNEORy7?gMwjll%hd2ro#mHp`<}w;p}Sd~>`)G#!*Ua!sl-Q<>R%W|`t*zq%*7vmJ2n*l~US^(esBdDu(YxXp zCQM!4b(iSu2WlnvlflF%CEqLBD7%WCEf?v3c5un2`K-h7vVBPXTjJ0&1Kzn~ke?yy z-`5rXSvjmRxD2n&l|sO+pN6W~%dPfB<}2z6!UnBv{Opq#(l@gn$t8xK=mU8_`XxNt z((tE*VGw`Kvf%4zMVmad5P~erRos2sI0$}NLrX~`O9pNR3@v?)l#I}~l33lT*Kc1e zMPt&)g!~MNwWOiN4KVo&-;^9>hIg?-qJV7j7IP*xX&~ukd)vd@nUmF(q2M{hj~(#P z^1)BfDnpLsHWC8)4PjZPU-*YFM0huFR-UENPj2E1Fc7uk2a#7m=wc@N~O#q;N%KQ9;sw4y?BT>5n> zg|wI(YBMYcVK%1c?%u?lB<`O2Qo3}f`hgR!p8V3iW7U> zyj$q}WeGE&{-McZOw^o5!=BO|nw@(JG>SIy&%IL|4@6IE&v#S%(EZz0OMNYYCSVLg z+QLC2A@Ho+(PZ_`1%rdGG80!FKWyxoZZU0si7*J3+?yLq-mpAFwHWX>BvJP zvGjX*B}Rou(8n2@iITP->!t^N^4Vj%2m4XH5T%iklD%fTSFxvnX^VDLttwO%ig;=` zn^AGvpZeQA9pGzSR=Pt;L}~wgB`ju%a3~e8K61})Eg&@yMvdq#Np7TfpvB0GhTcA; zauK2byLNxH+Ou=mY$ntriyO8(()m%BG6u2^DX@T{Vk%ITxgI%j_f6 zg7k4wfJocqwL-dGI8wztvf7J~;b3v|tX(+5kdWrBvdIrV0WwzTn7}g(K9Ke}V)GU*N5cvwU53r3nbTZq$ zHBcH<#^7! zYULXOKMZId7zxM%$nhxqNIuW>>6_w!reDUVtP8B755$Qu41a#JUCa1rs8wBEpvl)x zD*&^mNzH$J2-KfB$$&dcZN5)@BHPpzME`LEbH6YU%yZFCg9Qb0q|@ze6+77- zbNwwDiFT7r%w1N_<^yPe&?Hl$Ix~(R8#h_Vmqwg`$cB24>$C!7hSSyhVb1F$1Ys48 zIHSFN4-yh7kF7Gv{+s!h&fI|_)mOq5i={Qh>{*Gf)ELe=jD-T_=g%GoeL;Byz2;^S zky-S!;wdEzc5NYPt~8x6d3oY53)lJRPuxJv`y<#V=k4MRL}nkWMKnL+MJf8dw8Cbp zX>mQ5A$k-lPvELh7e-nt7&P!(o{yVGEm}@qDw&FH2i% zWezK)MWnzWEiW}O%|{@g{Ia6Cx6ifkpC1q`pq7zW**QrwABY?kLYx;k=xBMQyk3ZB zUZE((C!jF1*}b-9bQI6htL@ZtRO{!6Pgut7<-PbXA5bKnsWS=rRbOyBzjN0^O5HWp zV4be_&%pYBoM%+0Ny6ZKhz?yk;^yj+C$K`I{N68%yfh~4bgew(3;7a=QDSTy`C1U5 z4RlF7U9OVy>-R1O zNbNHkr|QbY=vNb!4Q$8GyR{DL2hB3B#sBmH1b_r6LhszYx4LO<`EWQxeD6bX8P|?Qs zXqf9zmc>hTnb0Br?@=;JIyle6Xg6$X_;j&;y?01qYNbZR=hdh`;aD6F*8I7iQ(_{QmAM`k$HP0tS_~0TuZTAb3#Om}(?Cj3vFIvnq)9 zfiyKi3bK`lTIi&Ma2VEaN6?A zmfmi3s3JQj0A`#VHC*7n2EX{;XjTbd(17RUg2BN8-#Mq^<+Uh}nK&Uo`qAHxBJ>@+ z9BzJgwa$%~5Sh)&V z+dhgRAtQ+9Emw`IPH>)b)-;T5SLCYU{Zc?9;3!{mwj?YciVwJci|yTDfFJt<;~1s+ zfbyMWdZK)1iGIcrgo_RhHgUSYdbeoDggalFWlNF07V0~5qt%mPK| z@G@-Y_L@BT_Q7F)o8zwP7aYm4yN6F?&Xu3ZgD~__PfS!@Yd>#_?DA=u*>Pz;8`7RM zk|#ea#yXC+d}6_p9%}Ud%cPq$)&aX|$OTl8Zn575lxR7)3%|>HBTM9}#Xa%c> zaP0SM4w>BMu(pr^7CLtPfLZo2=8FbbbR;_bys#-NkF)MI$|YTB24mO-T~aHV7GpXG zR`7Uy;14(|)J7o4Iqg9H^Dj;dNFweL&m4WFFz;t=&OTxg3?y92DE&mw*N4D(cc-H6 z->fbMa9}B!HS^0a(S)8k^j$*p5xcpQ)T+hk3gh%ux%wq9}XIwW{?1zo}qWwl!Oe)cKC;JoAx+rS?AZ;a)h42W< zDdA{VAsAG+>mnZ6X~o0P3n!!>&2GXT)uJ;BU@SjmE#3x|#PwC+7SEO)hIP&U)9;8u+X#$S%_hZRLSC$_h6lVxJ*@t+A2Qk|@5P@urLi zcKvVQ9ZZY8(pfLORPW*PzeRne$c%BoX-}HAw{}(O$45_F_(XO~QU7jjX>MMPY-_%m zVv`k`&Wh;We`mW0?cm!e(jz9_0j`=%ITFH3(es1e1Tc4qjD93K=(*?Fx*hNBusa-$ z>84S$#6}_tK*G%S^HEmHX)&T^It`dD^_U~l$U6qJW4T@%UL@bBGIP-;KHp8_=2MU~+ zna#ecN#{w!Lg36QcuDP^PDsn4qoe9@dLUFtT|eChtS@zwql*I z15+IY=r=k$Os@e9Ge}?cH2oxU3F)BGR})FTTnaC_9lxU10ji;KWlV27YXe$ST2bB6 zx2S|wmu^bpyToM*c=*%6X@0EZ8E&)K`&E~fW5mIJ$6I3)_(g!*^D-3J#AM3SWTtse zl%=9g;gT0)OUSF<9ngAEdDk8uTB%NW#n4uXY6jz{`N_`w(c5vRA?HfPbF1JOmusn} zm(FMfApA{Bvzwi0HC&hDcWLI6qRXd=Z8jUN_`pTYXt$QiXWJebIVY`eeyRHj`jR5i z`y#lHx6IxZrj5M0;67NTNT^2lIAV3l3JZS2fHv_2mh8I!^X&(Db2rU>@*jfQoHTjq zQPowl!{g|zz^y2CxsW>rony?>S0aXa-zXIrdu|!9ZlyYR>Cu6FeI1c__!3e0J@9P-=01dZHvWQQ z9w|gYY&yP*z5~3PJ7a;Jco}`J7pxcbeEiO-^=HD#D-G2JPiq8lcZ6*Qx%R5`oBpII zk!wu%RXG+E2iHO2jI=d=&Ai*|Wl7L*xyZEMkL}=b-oN5<_Uy9uM7~e&eJFfGPnVSD zb|>f3Xxp5(qY!H~#@cIT-XO_OUPP+?b03m0=o4;~E5yo`UwVKG^@pJtO-)sbBYWa!*uUGbQhgh3s6 zHV4dh+o#-9okB`I4ogvGl+*yw!ae!F9AR=fK>bseezjyTqrg$|mA@s#tSRVtO&6~` zWd}~krkgDi*oeSrKz0=L!#n!hM)IBfr4B`5(W{VwEUc+U8W4syY6o_C)UV4qc^;;3 zUp%r`?G(hg$j=4QDLr-33cA`|W?+9sd{?s|qXo0(yG`k;w7vDpE;|6~#xU~>-3o_0+$>-3e?|Ljo7z7h@B3~dm z(%WqgRZs4PP6Niacz}NJk9PgHjz*@}W5Kc#iZ+9_HMI9Z_`FDTg-gOlyr|b%Txu98 zJ1D?Q>?_J}0p_yZpW7;mn&Hi=xn*uXa?Yiw`2-2~1tDAXD3T~+#gOdtLc#e02n?g{ zdH{q(@0or|`2Ek{4|7ez(`_#*OHqdalZa}^oCPMe5e9v7Sq=Sw=`fRKGlaIzKn27* z9fwRKODaeHN4K4`iph1#(w%-1mcl)YmSKROkD&<(>g7A^W;R-!o$}238u#rNsWVxi z8!r?GX@mJ-elBC0hi1~5J4!+3iFHi{upOYyVppepScZi|u^zEGec6rF`bHusy^DG_ z{Jl&Ae@E8jE$)F!hkfrcchdjv0^Ax$jijoUL+viQsL{bd1e6lc2}=!f?u-^Vl_SNr zN$r7RC}>;2;wGhCT%V1sxmB7|!vI<>c=)E0>)cHV3ywKqjchYJxk4tE=2izN9qP+w zc4li!oj8r>-9RI|?plsqB~bsNiJlJHF7sV+M{rt{)MtZfPve)PFhR_gToxk>u6_pw zKHAIm4qV4*jKgS;%ABF$y>)erP`$2cZ4v&P#z~f-55ow)MhPBgUy%jCKRDBUL#>{h zbhu0-QHt8Rr?x^9N>ZZIP`p0&ij2aKLwmc(q3ZhTvhr^$Hvr@0f1+*+83X!wxU(0w zv-@YDM-4u=M4}S>%P*MN%vlA_#Em$xv%)IY|6rgJ0Pm={J^|(yhJ>QJ4os#@TovFf z4PNCDWFKPxM|5j#C&v6&1Hg34+<^eHrU>X76*J#>uhpk1sL1ObmS>{@cSB|F#$)k8 zunAc`*Y}ZQALg{&CZG1?^YszbOYd9-CVnCWzC`~#i{V=(ja9~6x*uVUvU~QguX-D| zg%*sD32Ftc#F`SaHfb;y(+NSMM+Ca$Y<;w!_6Y_6MijVTMlW^g!GSUJ2mxOVu-_2H z)oQQp@qxTST%Y*a^oore$rU!xcFyq!@0nHTl_PDJ*aSdqC@Kjw--3|y%Qp)nPp_qv z|4^c(gF_|-?+r*5ajv3FOg73KT+aTlwJ(o{^8MePEJebg$Zjf;HA!}3sYF7E>|2tZ zl%0_!B(xBvP?022_H2o45oIge*kv92GKM+tXZquO&+DA?`~7jw=kw}iJaa$y_TH}R zx@R6+t5+Tb0OC@3lJB$7i&Y_Y$LY z3hpMnjH(%Dk#F}fkqyzTGPpww+)|EWk%PVadsaMA^bnH7P_|uFE1so~?`Y+G8T$dE zjrqG@%we+v3l0IbaNZVyR^Qym*D_BntR`eYcWy}=w}*3Afe)E^=-}nb=GB+=`tCQl z>;}HnB1r*t*WT^vRSRF7pjj3brS1fy+gHB*p1H0ad{{q>>7^Nwp@0uRn^PIv`ZSLt zE-PV!xeiiThi%tpY;yO&ji=O2F|#_@kvwNBf7ozdqiKT9530 zxF=>h!Le(gLRp8I4g6la_VQ^a?|%S~ z6A=hA5KO%H$ui4Uud^U|Q%Has?n?O{J~cU?ct)h-5zRGJ#?i!3gOe&6L$XC>cbc0} zs%c3Xb)$jVElJ6#k5Q^Q%6l@|Y`uTO=T3}5F(~0TtVYtx@=n5{7^~aQoLY&esBbUT zL?1v6(=YBhE|jBKadLK0t8Vws0Wi6~GwCtR`;_n=CmvsWx9(C;H(DbO)}z!wJXX)k z2LzOtKicHbHcEg-H$~MPD`+HN_d2Jm-FD>8T5}NTL!kh?6lL~APoHyNna0r~3sn7) z)ZmbV;zX&^Ys-iDUQ(XDm~`S0O)4vn5xFMNPRvv&zu)-wNMvKunCz?*FLzlNzQ$4o`n5#yll%`)r!(!8y5BlL$cRm0jj%)XJ|WZYpY zl+vEhM}BE{_A3-NFWWOi=SqS$n|Rw-k4C?N1^uX~e;4SS8tork`uzk;cWrnZnOdmt zPr42(v#F(l-Sd>GF1*g0ojomlXmbjKSbaprq|0!JhBC0icM0L${ghVUOa#KSuI>I? zT@J06Xy?7> zL-zfPPegHAbpF|Ne3?-w<*g4F~d?MQ+YZunLv3}G4N=&046P}l1C>Q==iP4MC z(a0}1QB|J`4VlP7$dpuR^xTFe2ZY`n07c|ifM zbe!PIE$dl&OazYeE)0TN87XuOV&&j|u0@$cN{m77Gn^pe2^s9#P1%!$GuW3m{r|0- zzpwXDLL(!r!r^|>3orkYz{{Ihh(-X!*=hgU$TRmW3?4vkzjleSt$=0c_TzKAr;+7^ z+L_XD&V1sa45RP-R1LYW$#hQFxA+eEBBjOL-mwoux}X2;Mh4%Nv2fSH9on}BrM@Dt z;*7chddie#?LDgu{hh4e;gpeqZm4c_kZ#=nBSVSQajU~~Kg8aRA-%uX&}xmSx821Z zL9A(8R8MJf>y@2W^}Bdv$A^EM_REv%-MXTs(!J_IK2Dj%^Clhk9Fc!7{pqRjpO z!>HyY(NnB|J1!6@>pZ8l!cA$zWSmiNIsFiE&Se)WrG1jC6o2`|xVprq)%ikl9Dk4E z^meA|$ejMl*)AT8KTFxP8P%b9?lNkm*Ng+|qZG>>wDn6{CQxQqzhk7E8xZ)HV*CLTAEmkHWZh6C<%YL3Uiow!G2AGi+I5%0s|@xavqigIj3e)(Oz_nWox$A` z|0XET2St7bCJ6b|%51ckk;K9YHwvs^a6TT8fZznk#4RI(5w`1?&DGQECu!ijw1RIr z$R;N(u63;6gbo*JXqi{9J1!D5_kvZMuAuce;HrL@E=$ViAM`McD!qXhO5Vw6Y6NL; zL1-b{=Xzmzc8KID6qKGAush@weKLkpRZ)VToUtoD%JYL0lv%ahC!hS+&XR%h^UWSJ z0t@2|I3ReP;lY^4U_+=PTx2q({V)D<$L$3qPw7)tT^RJLsHMC=KD4JyObikVG&Orp zX?fF$_g9O`mNaFC+Tfxn_#jb`@ zC;#12PCsDUsHat4 z8-=dU{lt$G4wUu%FyH{c%AlbMru=$vtioiNxev!qnY`gcwD`YP_X%Ybk5Tw*PDODS z$ISAJcL2q~+}%w^u#{I1YW_a&bNaWtd4m6Db{hR}W-rIe1Sv+)b`zld+#mMqtbtSM z2vP>Y=cfIX%qiCu)X&^vK96j!uneUO!7q?sz+aSHZ}IVa*Ls2fItE_)r1Og^O8{0Vnjs;JX&&GtgOwmUc#;6~Hx=@IxU0y8m_{ zV^wn0y9X;+IZ}TL>2mEh4Up5;Uf~(WQ+zTcQ~~W&9Hqrx-?lfpqvkN6Rv% zs)0a4C0?bpyD9v?wJT-yLzXM?)|M7YG4LAPBIW>aW0Y#le@+~_LMHUTjB6kYIS-}p zEQrirJCCwg|FPO?!p7>Vkw_AQM47;D6SsTWOM{&I{MextC3AvR4!qIwOvWXQzxTZ< zo6@7Nl4+BHcaqwJnjZ;Hr3Ehd)@%bN1{0_Jk^TLCcMKFpGw&`SOw1C9`c`(X<#gO_;&Bm)C=Ds z=vO{)-F4UVNn_q{;_G+?ZMvl$0uDP9`dzZdg6`+w;9zgB<5F;V4PmA`H=N%j-|9u5 zwNl##Rb4x&dV7KUq)jM#lpDX<&{yo@z3WoF)?;qH)nXEXuDi<^xcXCbtHPXtBZ%h1 z)=D~=XQ1q2vVOs&Ur3QmbAp5<2QcB9B3?)T3)i6Ovrn!GB*y>8fmv!`>a>Jaw(<0cNZrF7=c{%X`y8 z^(33MnQrQROgwU)HY=P%K`U)AR$x19{H@MQzWr{cCk)@5wzj?KiEO<|^>UKeSOX*X z8P~<(TKk6kTg%Bcqz(<{6X&w3tE;yOHRR04x?IWutHZV7QB1rF*IQEzb>$r07OiM(rh@BVUguD1BOb!LJgu`6ku9Z*vK@B4B&Ilx>#g;uEnywo$xcE$dAWkA~Z= z8JoMU{kkYwTzReSnGE6^+*AAtxBkMeaAW!W#~(6qNF-hkGqDa{>(IbIy?h32W?Oex zCPZnaReE+U;4Q0|>8J6d{96~sn&7%yakp_B$6^DwM-fgNhLWF5vFD0%E1!7M=CRS8 zAF?)Iif5{tYp}>yT?xY_aLc;Gws?SF?rszKl}_F0$Y`LCI>5Qyh~51;Dl>p1@6-z0 z4u83GR5RYS@>+WEx|42<(sI*L@EW*xAVtcatCthae`I!zsH-@CkV|Vmgh@iBJ9IX` z#)?kL_A8Am2sDSS{3nuZE({%>6uVkr(~+Up?Y#EKG5v}{cNi^07+v)@c5{o{-B-ZT zEUc}qtE|!$eHR3iEm-UShBGM$U#35e-q2w-A~NjPelCoYUE9=GknWhNm~A1E8q@+kZct&E#0of z_I#>*Dy@PzlN;U?GLw^CB3kQ`Usba5Ly-Y<2E0);zw%TW3cryoZoUq29t5`IP`}Ts zxc-g9+BKx&TrB(i=h+YH(MSsGg)B7sEf!)C{pD)6JgxsRILzMC20OlFmSUfccRR+XFK{G<%@tL*G-snaev$@Z17 z+?j<1ro3fw%g?I?bxbZbuI(!U0Ra(ch=1Hg_dKb;yI=OH2t>%u71VKYboBOe)IVen zSq%56S~K3xeIjnop4p@)KaseeQViSJJ-aqJ)|+RJgtVkz2@E56Zp?p(MUw{;AvgQJ z@Y#yZ@!kS3QmLMl-*i@T*cPpcm4DVt)}S$NJX)R9U$bMbbGM~s%4Wk=T<7U2!I@G3 z;Mou%*TD3}X7gK9>FKf_R(Zz#kkPCI=Vh`U4?vo$T~sR#NeF7-(8} z{)fnFB>z^EW;{g6`SOW7`}ol6u3+l8)3+33*~PiexLz$ba|?PD_q}@JiFscFywxA> zNPxn#?no$~;=|Mr1ri4R>KC~_R1;@9l4&}xGbyj?YY5V*FG7msSn(}gmjA76&QLlO zkc%GuLd|@v5){kPC@kqEdX&L>^Q$9v+Z0nTKlS`|z+Hhm-apTAeB<}|cF5$INNtQz z+Tr`k{2`ko9Xc-X1yc*}OQe3wfaxst=&^W3iJ7Xf@71ICL`U1b+hlr+EPNJWODgrJ z4fh{^@;q>%{M6Vec{720qkQgqdY5ZTn4`wICD;X~WfR@}+KnA~?@1a5E10P>^UK}V z1&R?!PD2L*YB;Z#=qGfYiD^rc>Q%NRPPUUmMk>u?f za_;I0+KI;bl79DFVztdRi;FdATEaV_@9Wpdjqvf## ztcpv|`}UCaLQ&)Qa!?_%FQxX}{N1UCeet_ShwdX_8BJcaPJsV>aqnpRN|&Ktu@dV+ z>9CjSO0}xnJ?0)sTTG;-1l+BqD4$!5V5YbZr(Z0hoA%-sWd$AhH^Z^$RxN6Ujzcj} zLQ6=mBY(cuY-c{2@`|nXOZn`Fic$Q89*+3<*Oe;r4xwiHO{zrorXA*NUGZr@Ku%-U zA%v=p{TQyHS1+KVH(B^>bJNxk27U2)tFt_Q!*^;ENt<;5;L5Yp0%p?5r0iJiExr$! zyH)Fb7wMLM+b0|%_OEV^CG0Dfpkm4`M4dI#U|3$ic4|d1B^ym2yyzYN!wHDY~x%; zKBY1##U6LQx^FV|1K4+m#xmNPA=<9~o$GK^S>eRgcjXwyh;-4nWYNTJtr3C=ff_=_ z)qv3ACZ0)IH{o~JOKY{$Crve^9e=o{74=r84Qucd^>2I%`T6Qx?eIn@HB+>y9(XIY z`Qp$t?MOUjYA@#I}v2KI=a1L3snY(l!OX&*v*OdV-2Ue)+M z*p4NW%Dl!R&;XFrUeo;@NZ=?AKvg{s!y#$K!?Ts=5iBZ9CwaAC7(%@U zJ*tS`(*1r>Wi3aK9#b;t(+_18DVlS$vrwFONn3{5x{R%jKNJZ(odv?f5cG7$-AtUP z#z$i%(-eK(tE>pKvw@J?KLGr*5v)zFUd>u(b{zQ{S*tw=0B2`48yh_N)G;i&2I3kB z2j;51l5^Kv|CLFZnn8f<;aScacMp_9>Y8d#=NEfpS?RX~H(Qs(07&8*f<^ z9x9x(3bz#8{BFJbtGipqZ*h1k9artCeG+hvTKL{8B1ym&JR%MzW@Kuk7a>A`7g+b( zCn(=#=2x-#63Xj4uvD17~w@8~yc+#Y7s?hzC%5m~8 zfy~7y1K{rzAoa(3DqECD?gjTEKo#5km0W$WNT2ugjgAWZ&H|gR9A@TI>Ql)1usRI9`1~ z3(vSmCamd=^UR8Y6|#R!_i@kJ-0S@_qtFzvQ^#mxYNnxSh0i=tD<$vAR#@CRq*(Wb z8>!~TA1t=E+Ht$!EPa4g0ngXk)XCVepFsU>ciJ5!gp>lR*xB8vNtg4o+YX4Ylydy> z9AO&B4c_g_9u@Ovu7OgEVz_@;YIG!uf0s0cC9ihs>Ag^pnKD8A@JFgF!7!cJb)QrE zS*aUJqlroF*tIC_Z>>+It)X(-WOonT?7Sp94;8%j{xvvZgWljt$S5&*@B>hJ_=o%k z57`FQ;Qz-3fncR}#lA=lc!gkQ3o<&T)g=Z}->8ZgO@KHX5 z1iv!Ni|3MgmzyiPxsiY)*{kr(lI9$-2=0eu(r}-X4HuOUoNIZ!xmJpL?}6VyYuw!7H^d3?$&?~LcLYN7xwU-0VB;m$O`>!1q+;jOCz)&M~up8Da# zhbaD#l!PXhT$B9xnmRy*z-yyfjv}nN*rk_PU30zd8LZN`bkK)QeGc#xq{ z+Y1#53BtMeybfT=TkVl({>?A&p@z1c&DQ2N#?EX#>{30OFARhvp;V9Ra<-JZre)QW9>I`_ZkQg|cp*R=_GytxaSkb?PqK~QU;POE(_yYUJ%6lFyFBqd5Me|> znFGNIHc_T}9rBPprXY{H=NYG}1QVf9Y^oYe0;6P8C;ipL_SwaXo24OKLEu-fHBOn| zXi^Nw*8>eEz>b>hzz|?b(apE7MJYtdxb!^t{7OIDleay>jFSV(t@2>iA1Y*hL-X`J zO+tG*wkAJ6pi8e5j=YkNOWLe;8y*6HZjrI{i|tHwB-!OTMkVC$-_gu8+p82#Y+cH`7PS0uq9@^ zzM+#A2nE-k31HmR=6fEg5K!^>0$m|rsA*+Ha!OyP!$7_8O%xEl-kuf=0$0^i12-7% zD})ipdlIVR0dGr+*5l(^J)M*7>FtdZVQhLzm+8EI$WIluvs(h32ZhWHQxBo*uku#M zZ?SeCv;6o2*&)9%(7r06GvqV1U+A-+-$bY_huGDVUZeE@k2dZQKUh|*gP~L9yZ~-g z%fZxas-!Yo0VnzcKFv2Po*}KYlc)7?+V=Fm1r2h&@Q?Qid;;5{CDGiY>I&V@?(Sk> zY>2j|7CF3F$6Ny<)%iHCZ`p~vyegi76 zDbZ|8diL-C7dtMaP)Xp?#UQIgMS-X5!@&gTw*TdY3CvB+?WQ*QyVWc0!HqF31#wtc zJ}**_9sb>0ooL1OyM^?ZLl`y1z!fa;@*nLLs_$x)U=e%H{16#M6~d%n0MRv=#wHmnr~vx6GTgHD4UP#%AC>` zIJWSoJ`n#`sr>WdxRrY0+7QC_%KJYt+7Ax$nB?%hCvF;Wb3FdKCR{~0FX)Glz_y1R zxlC{55%Jz#HCy(5Uaafou#Ndr!kj6tr`!v!XPFKZT$Cn05ff4b6aM$66jqC;H?=w7 z6ta7X;Xa1#jM&M&{|^wO;QY6ZU=>;vobn?+B(D4Ns}MCxpGo8pasmG2#Qzz)2p@ef zwHF0XH3kPn124gCeYZw33J8TfCvCn|46S{zkC&f>QaCd^J7!OXIynxZmA@)9a{<4h zdH?X`m|;kR$n1@y(@FpGIC^ zSI{>}m_4YXKhC*%piBxy{GFXm1)--Ekxwmv0erbb?c3*+z3v#r$bFprrt{6EZDMat zH7LK%#h-zt>>z~!I8+FbNef^-j_JWZ*d+vYLb4T)`0bEv34#6rs`qwv4P<%6g))1{ zjfd2EJiqKoo7zI5Y+B(sfqO(Y5%QiU-O|s5kgmlz1ny9<&zBs-><^}e$Ngpzaqvs% z9^lrX;CYT+7}?$h4O%?nu#9upt8=ktxuak zPF7)|447Q9rH|?NF%IjXkLXc={ccLX;Q6DNdQQjZA3+2uk&6X-vJOBNFcC;jQ)SJ^(C=D0&#Q;$v(H7QZrUGMiORT`Su-UT3XB;v4?-Pyn8prfT` zgw?jdX;$xhi}IU5{6HkJhU5)|Dx^bg-#kc-nt(vhaaa~zYi7+dgeCtcmMIR1Ki?oC ze6W=5t+wYedKDBZ`*xTnpZ{4p7KO%JzYA;ueI(;v1OIS^bgLfS{8$J9Be}+T&p<16K)`Gs=;NRaHY3n3A!z(4fB_|7 zNA6`Kdz%dr5h?*RJP?!u$KqslQ%>4~CSL(U9WZ(Daeq*&dBEQZ*P=ED_MnzEs;uzW z9$x~fS_u>wy`WJxs23o&cmzYmh9jyM2;GiE1IIHhu({Kk10cbx`-%x@s&boBSM!cX zTm>AmI^7A{bK8SC{rlhG@nWb4ePX|$uD!O(FU1b^_HMjn4!%I=wf@IZsQNolU1r6( zlbCwYscwQioo`Y8Yhwpg!P$ih0uy`LZC_ugXYXW?oZ@c-ZvEakcRs;HBni}SL%Jn^ zIGCklvICxV00}_3RcuXD-nQr8{ta9$5)JyZ3Y;A9;YdrCK7Wuin_9b5n6oCYB|bvP zbacvyG+0_T?M10RpH$R7Pc!3_(#1en9o4fv8?kcdY3;X5(@MY)jB2pPK z8h7b1P^9BJI%&%@fq%F{7?9u~b6^QM_W{VsaIfJdjrSc;(OHv*pE{0-9;}-;ugx42X3?>6{WL0+sbw9tF=zMX_e`|Sy5VBDbi3V1} zeFz6`>Kz3H4e(V2n!GGYX3Y$JR=zcv8j0?9z71y+Z3+3-HM^8^!ABRL@gI1?p1>vJ z{O6N{_>}2KyK@}`BE;X&Nwuzxoq(zkIIY`{0Rtq1<0)Qso%0{kU{)>yPD_sgJXPpE zsIVJ?LVa+2r~+Ca1Q>@p&&?(bLUR7;C;-I5O*Pb`kz7M_+c+Ex+L@HsTkJR&- ze(>_adRZ*>_fluN%Sf(^J)ZIS22-w1>P$(e5wWC(LNOGXycVvhj6xl7N(~#{ZmcLq zlqMYU(>s#e*|t=Z(C4V)*khUb`gP8WtvSDftY^Eo&_1Py2<2a~tKcK7P8uKs_o(ys z+JsIOsyM#cl$q6}9F=Hx3?;`YR$H^uGW_9m=(7`tG*iL^;uWJL&zcVfXi?`P?>Bg^ zjD$TCQR_Ah3fn1uMc*|9M2>67MP|G3}a zz0ek>j(uGq9m=1YQy|}TqB7-%+s&rmFR$0_ zy7H2+*sSd*6iX)a(6lu?TBJsLb*#YTeA80_{usOH8?J_q^#?OsF&Khla|su%o;@_T zOwDJdd6*fa81E&z@|88zb4ICzMr~x4eKA*7Lk9#>tgJpdq z>Oixi`(pQSuD$IqZy&IxU3Y4ElP<9`tlhxTv9|h-O)Izj7SRC+cHfTCF0x4H@V=K7 zsy>5BsAExa9*qk4x<^fs9G?Gt&0Tu!ZXsT_k(>HrtpS@N`Nj851JSN>_z_!3)M?7p zWx+1PNUJWxm}a-t@%41|}~53x}k zSg$6l?6cf9B?==C_h`1n_aBfQeA0QEvXCBAHx?&8r&vLwn#~Ol-)9anUbp(FIIf>I z(I^labY_fkJ&eF_KpqGLN_S za`=|?AgJcfRmYF#8#lA_CM;hXDj12d_P($Zgze+->7TrgzgyH2UFS6Vd`d!kPO)&l zHK4B4XD(E`*G_8rOwnz=E~DP%TVrtsK?`qK#eOsU`p0m)ib0Um^49=|6R(erl=fC7 zJ_MJwoh`{v734FY%+<0jZ2U5m;9=$1p&EeWZWbZBf~nQw47=W)9%L;fl|?ROw-fjTVL=-oihF6eq1(^FVZk# zZlS|S@4_aO@70Pyf9>u2&pa;cwfjEHDs#r2oE;jfX64oI6rOi~F=MKdL}bp6u~#eY z7aLifjQQ%c)E^)eU_SeFft&jM`-SYVVFRL-R{@QhdCk!I+a7#0fy0NKZ4gwn*wgYA zLv#P)Zd;@}BymFT#SN!`u5j=x6K$XQn9=W}F(MM%(o|K>dV4rFp))wU$Ca&8R4h1$kHm9zAQ#e}L(-k0|wmQVkUPtmIeoylrd^-BQ zNvMp=BDW;E4lg7Ay~=}I(xq;vkuF`naqoe5lg604DIel`p9C)9Lqzd{+B+CNCYiE;H6srAZQsDk+a_B9dkAO~rdwA^vagXLs7aljxlZe>Li7B6HWZ zdFKC-qI6sVU9)DMW!Vv88l*~Ys;k>v5$dI1esW|_f?`SYiN@bAH{!AJAL8R22|8%4 zw1u4!+_ZC1qI~8`etNj}|AiJ|_%w(4E! zLC8Mj`M~2!Oee<-&d)5d+-*E{^e!pJpjwmkZ0n9+?YDr_Z0?}k&1W;7qOBL<)tvD% zja@bISp=>}d=HgdrB4&gBhNk#lFFkEyyPWR665+D6q<$`TbM#s2}TT2>aZ zW`Y@W$BXtTOc7*=WO^7{+chjHqkNSs_adpn7b zT;T({n>~lvnW+bH=_GF@y7i8}B3x_vx_@bFW*$xqZC;BtuUQY_HIH)!!nOsM+nu_< z@OGO#;1f|GNwu>UTcm`)8%y9-q7rG=GHN88EH1y598q@0kJ}>yxKAFRw&!+}xQ(}C zshPI@;S=4K>7(D7Sa)Do`Qn1wD9U>eIC|B7Q!`Ng$-N`LE17puV+`>?5B)A zy-09!A9A19+<$z6f)nx1bSTN^mH7HFEIz~ik*KSkX&>R5iqC3Cy(n4C*(&QItE%(L zNTV!e!H|l?UabTxFNe{n^_e=&{lbdXOSpPA+B`OG!g={x$sgMqBdb3hsGcqj8xXIV z($9n9#eLd3G~dU$aWcuG-}~Fe%-Cce0R5c@?>tXtQAq0J-bRh6re=}!4>4VzU)4)vFkimRrnt9P4)YO_kexDp@h*m5>Ng-ZZEKbC&OQE4bP5C+&) zOX!fCbUgYz?M>GjL$#($4QQkUFUKY$&s2}O#^D4i+I)6}XG8?6L|sV;J$PC(P=bG@ zvuxF`M*tJmUKDAb;H|lzd*d~v0I;ty3)czB&QmUFm1N2fg2DZ+#>%?7b0Z?2GOEQd zMZE3yY(v$yYa{iId1NYwvv!u2+CjIf?0mHn%FF{isa1-`#c57m8~$zkA>Cb-l$rXK zW9)XH%asOMq>(1un$nr@9;K*~Qh`Y18)X;tm|&EQ}nL96R*`5@xe>&(7?g z&vWHNN`)7D`kKun`O@P$&*_fn$Pt(2hfWr#h^UrKWR(~gkFeD!s$pbSmayt2En-G;d&s0e zl8m9U{k&{ce5~7*INVjEHdE#FA&aNBA&cv0hr^m+jAAI>9X}X(Q_7*B$2eO?T><;U0Bz zF*X@~+DI70@t6;vzMHoVd3o{OR^-LETZJDWk3aILc~Ezbwi;k%3b_YgiRkfo!c);I z83EVlXViR?XsQuBV;fKoqCMxjZ&Nt^s1Iy?-XtQ$;M>M6L%)3cLd=p8MERbxG>#Q^ z!27qH0Oozq>0fscz4PjG-)}dazOkVWt4S6J;~twM^UK9P*N;U(%HAa~_zPq)-*X!n zSK#=t8X^PZ`ectBn2B*6`~}en{rcx{cnJnFuJF>6J;LB{orlgI&=XMKdG%M~F0)4)HdhY&pd4Ik086cqVzq{Klh&(-S zdNT)Gt;hJEo=Z&{2r*C;zgyS>^>j{e{I&v!Z*sAO--3?fhVIBW#;%f%u- zfPll-TrE-urIB!Md__DAF3hiTx^Mc0WGO4q;ympR>Q0)y`7Ok+YBTttZbe25>nn$J>Y)--0w z)}Ebva&SCdrBzmHWnd`;vz}SNHHx>lH}fWEH&zl2p(G_}dfVxsQ{5_o3=;aUkR}1jWq0hNt~&S^coR@2rQ7D8v&UTO0l=#)%C6!u=0VQ zb(5(QT}B3wtaYv#mY5AxtzK^GP>;OfG5~D8{-!@Ge z6FCO%BvXlOf&eL#8wQrCi?-)nT_a&JgNrAwj0s@YEr%qkCk^=;I3uU3DYXd_$y>t* zIt|N)XB#Ng=Xk7aw_RljT1&pB)sy>rmsQ;Mec-@yA$fU%lQ~+2O@id>dAX&^`%eh3 zbYWHb)w$MAryC|Yu8ahPejCHvQ$Iu+=M^pu@VWL~g^GV#YDW&iA`^n>+pQlrv|b%| zubHCh30zaC^l3GZ49}@@oTD!34!LmBmavgFuPx6`tb~WeVtpFVGm{g#qz0}!*6j=k z9q28X%JqTO2aYLHP{q4NImiAcHz=+_EpJ@bh$8B#bI-K0zoqZkQj7%z1TzvgEin8^k9y z^+J073enxZP^mGuEO6@fi=P||zX?P01<7<)onCh{YIwFk{HAyPP64jGiZrvO{?NOxyi0lXQf_z6+J@G@Mk z841_BNBP1gBnp^BT~%&W5#htj#=y$%bLdxj&X#>(fs^MkY}aKazbp&o-HF>tNKw+sEC1b$GOa>g8yl7MrmhHbuP8R! zmWAgO)PKA6PC~G@JZR5YEb1dmKhv7X8jVT{$!gMEab0K|8v9;E?#y<4x34ip$ve`$ zGe?k{qBoS{le7_uS+zo8?MuWz9`SA8A11wYqcK*3@AU03=>0U-tf`njSWS)0Df;R` z>CKLHzD8?ESI?h&ClWX}KQQw2{_L*;7x|QZCGWZq#xIdNH(X6mIBlWUXl-m@J~tc! z>)|!yc(>w7SEa}uT<bG9q@OCSzxwxb&uW+u3Z+^^x z?I^j$o_}e;z;U|78&kVQAtjxMj2N0btlg>3<>R(?+;!+7*@Gjmbag@}?{(I0A$3gujWj>5EW~jjeB5(o z_12F|1||5i+znL5Q@m*n6E1Axj0f7ho=qAxapup7>CCl zh;a@Z4{JYcy0rLq58q^-K@GSS6RU8PST3pXc#T)1IyQQ$%IWMkyd(F%#G?*_vI3|p zd!4Jl4s93tF@g#HJu0 zZ3qyN!q-1}*|AD%j6C4JF{ivw5K^Sc1;NxElTOys9;@Sv`9&_qGj3gm*SXW2V-j)-o8J1&Agj)Tn)O_>)NuLz?0DH;&vG+@PVb)M4cthx zwNc$ibYVpzitE5hwb)&(wQk?vAI;@XkntSEED>qdG*j%Fsb}X%0`=3twEuRIE2W=r zeQWsGoxE|uea}I`s~ch_?afip_qpl^LdhJdf z!bJ4B4Sg4nCa$!G05e{Xo1TdTuHDV!l^A{z;?z2hSqFc$sagLx>WP=(Gq96bB`?Jg z!>Yz)8lQ<=Q}I&$oV*b!FiHZv!%nmPwFCI(fHD10@xH6i_B4KtRuWD|cBv95D0!GTpJCV;1^yT_5 z@+oix-*7iKz3B6M59F5Rq0q5-m5MLF-LyS(6c{<6r{KaAX76)aV*bPM@6R-IYkX)V z%<7${bWB`f_igmQHqfUHBmUN(FI6+hEm9M;y`BXqlUkaY(XNd@IxG`%P# zR<^IFvf^y3EZ)FIYF|EUgkwgSMS|+UY9~^XF>_xj8IAU?p?$Awbm`?2sd@n!)WTi` zwh^=|XL=6EPQB7jPmo<5`Ip{$Z~E~dnFyO|_d@6q;Fe1Jo;^hBR>-q|d@Nb|n4cV| zeAfEfT%-EpvoPMq-3lqlvi&Icc-bjq+?Wj7-a~8XSrOjuQp+74^U)7r1FdTYxhr3` z3Eea6v05jGLSu>i;dCS9eVwzvopaD*L!mgghC6l{Hm;0`SD#z>NF^qm6+-eZn-r#1Y47p7chOYL)s(TxD$ zC9))iuOZEHW$$F zx>8EgmOLt+UiEcnfJpVaZN4zzs#ecWIqI}P#p*;oTji8!ZNwO7%hxVk9{HRM`>K+t zY5n_%kzj>Gkvn7=o$SK9Hz5lh?`pFQ<$Z=PXoSXt^m?&wB(P=jum?f_xLk?Jc~iH_ z>FYtKt*<$`B?{T+;#{1T@mQ|zLJ0n*RDrlp%c~=iY?zLjR<#(|Nb2}}&iK{gk!b4> z?9n;;G^;J!X-5Y$1@$cT5s7-cSnmrxhWzFuzXt7wGAbrs1wbFn9TB_JNh zqw*x6?Mkl#>0-}Jg#r>4Q1xfDM?YAy17L#dYUht1Bc**p8*gK@FT;pSL?r zWO{RRHyH)Pi@iucmuqAs{2k)(=EWYQ<56ca5k0Oenc6pI8zJr@1%ne^~I?yPc$7=+B}K{^6ymCO|ye$A>cxZ2`qwN6IWf~HgXQp^6x zy@#OxpuiOcwQeI@P5JP$@R(&)TVbWjj4 z9hy^VM|1EP4d#up6)R5gUhK}CsSp^cn%6Ap*ztj{GeFDn$EaGH?bPa4yWmy_T%T!9Nm6j1C)%qUpM)^MF7L=O)#NqSh ze&-M@dRo3hjwGLZCNgfouc|Ej{r;kKU0bY>XrHA;s|`!d5a$n-y7V4T9s9M~bwI)d z=~L!hNf_{mjj^*@Tk>Vn#0_!;&3#KoTf~JIj>V&yIZiH^doFg@V?7fxh?$o5xRQdL zW|4%c0zW{x?{2GR2|>bj>(r?Ng*jtdIFyHz4%hUKjHnnZv&ftw*p{XC`FLg-`6+}~ zIEDV%c`TvrNK;cKKj1!7zmB0dQj~(jigs z(au+`IW|?fqsxf#ea^Z=79uxZ!ZrYS09xf-igijy#9w=hBII)gW51v>oN^Y{9;9*! z5{4Vmo-^IcgU?TI<|1LJM@MR?Lfv=k%*Nuh#QI)21U7dX2rLc6Sj`PnMP|B-LMvyK zGqZ-nv?A90?ve?a5yrBFiR|E7f+N$qhT zYis$=PQ26CNEN_jTBNKEvNevU3#7-kmny<*kj2%OTU|9_<8WoUK?}NxZ*^-0=STUv zZiGUAIBU8r-#A5j^-dwFy3s;5bmDXZ^3q3sTg<{5&s|8cjI#ifB=*_A_7??)n}h(C z(hC}n?Mijqc0$D7lW?&nF)isgz#^g}VfCb-(XL9zMUBQSARLzTxpM zVPwKP#u+>KPHx_^A9;3xH&mE;l(?~g#$irT!T?djz6zxh<& z8J1;=<4WP}PcgSpD3hqugEE50mr}yoNi-I1o=vEiwMZ{$2W9zJPp!n$N$rA~$|Ov= zrY4_*#sSpHwxgRelrnZ4LZeX4CMnth0l)L9e;WL@C$Z#aTrv-w3As*ETA_$t#bZ6~ z#k-AuYklMo%pc&@T$R0J{^r<@@sA?G(7k=GxIvOz5N4B0scr8(4(~l?J68HMasKrs zXr3zGci4?Q+=b@ZRIgcv5RlfiVSUq1;DeEuwtW)QhGB;_>rZGZqtUe0eDqlSB-^Qv zlB?7uEf$Qlc^i+Ass7e|Fw>K;Qv;z0uvvXxRMFOb4l*Ilh$f*>CYD#ZuZ15X++Aw= zEU&yWDjs?oe+Bl;VZ)wPmbynhyCQ9QJFZRUgkFh>nD}FRDBinuyX$Y+Az^B3R~K0lrR>&A&h7+$ zuCb4a6Jw{!P3JPzQ&(qPnYpoTze9IPbNzLJw+?l+OB?oPXtTFlloto2Dux1~4HtSe z=!)3DxhNT=nWq`&a^1Abdr|2`)kz)KrEm0^i#sA7uqZoISI}(6u~sA;?n&e$4a+g^ zRA>QTYC6n5291U)axSvyLG~Y7&N3&A4u@)SI$fEcWHY`}GU)H!;6B%A7Gn#i9pW4p z+A`H-D`;LUQo@8N3$e%Q=Z+9`t3_AGPCqoaZYUlpXb5Q$SBFSnu{3aG`q8bmalzD- zF=GSTcMjRjAS0jDM20pTPDofA{`fD|;R5&UKOL!9wXI1yVZff2+i8q}EUFE8+H0}M z=))du9mxeY%b;zjlSUiOg?UmTY!^S-O!Gt~-NiU0ccISGC%aU#ywMW^`c1>0o_fiFJ=7371*g zS9r>6GB>ldxAmFlZ{rj9QlLv8w!eC@L_Wm^M@|e5;gukN8cvHj?! z8;)wqOy|&ts16RFNdtwY8;^@@Y@{>{+nX`J>nqb%v#4 zsxiFX&Rz9bkxJo_@pRmKx?j~}Q|=9N^dG7GaSX!G@qK%eA3p#HVq`$RC(E-#ckp%9 z6tOZqvRHBH4V%ojmtCqPsimpe?ntf=xP z&oh6No+=zAfNZPP1xRR`J*PeZ3=II{5NQPY>;!S=VHX8qp=1Pb%3S8Inj zU2n|wIR)%#^<@7e-p{Ip09}LKRQ{z^y-yfO;YHLAG^er4?h6+Xu!`Cw|onBUA#4%}?+y)QqIGbJ$5Lge)Qn=3bj8Y-cWkV--3moE$x)@@(5R89sy{?U z@xQ!DjPl6{9&B9z62%{qO&ko3-B*y-)DY|Sr2J#_~y-* zV-c@E=jDOL;$?d0`{bSS-rr4+lPHp?a>MUx)?*l|xXrFIb>YTuORZ9V7&x3EMy0f? ztk*2XFP-tj{x}~{FKJ*ASI_srGr!?)2Lu*{4Jy-damCtgr8csCt;dR z<=yg`^fLh7{myj>=kX%}V=XHB#Tlff>V?#LLUqV6gh~%Fs&Pg2Xshk>${s6d2lnWt z4M6$q{4FW-M2$4;$I|DR)_3Jz)yOL8TJ$RrX85j3%5$k_OD5V&#DsQARmg_0NMp4q zagmY5UMNhhu2kN=={Qq?b{b0)mbJ<1xz}N~x)7!+rbMn~YJ@_!>qbZUnY$sPSRXjD zccyjx^SzGY(HJeVtZr%pZb_^@L0W=fS6O1--6LTeb?PX=dnxsvGoLZF5c7QNo!V@l z05Q!QdevWDKOb57rl7Sxgem=VTmCj|%WF-)x@ZuYd%T0o=?WYNwCnpkms7OnAS{&B z9w>u`O-`85#us`$g(DFG6hG+PWjS}{HaNz)-O9?;(qe%MTVm{!J+1}1?(D@gui?C7 zakaqqAtPM@dNr~mujg*v@Ir?c?)wGUSRg;~9ZOJGlR0ai91s$!o;T)CY0^C;-n53N z?L{{$GXNqy@?Y(`Q|0(y*)d9!#W+%mG}UcxHb}g+;XtyGUN~&yS?4L8sNf|m_I8PA z!lpa2DStXai5!JrW{UcfrKvoADdn63B52*`Eym?1C8l1PHHB$YNTf|4tvKVHug!hD z*r+5>vX(G)<&}n#rsOO?H28d7HZ@%e+M8H&!hRz-@dhhS!2u{%gExLU@0w0sRzU}H z{(_R{KyTa>VQUx{U4WED+cdHnv1EhQK1MkvmwVwdWJaD%89EHnQ6<~x@qBOMHKbCk z`ieX}5P*m$JS|tyg4%KWttX>c^Sw@o<)v*OmFi>{RDjIUr+a(1*WN{({{Le-q}T-Z z1%FU4K&nXJwqN|GHzn>ts>s8dzSkfBAS}=i#ywz{^7$G{GPjwXmHyDrUz`grgF0r9E1qoc=BbQT8mx^q9lftyU;(R4hKYTXg8C<}ce`30qEbA#%Zl|0IM| zFVnxM`-Ag)6Z3CQUxj!G@-i@XXAZO&HpmXKQcmBmAaDcE! z#{t6Ld00A`=s(??V2v7500hD88kLAk0~iq~eAs}D0w_c=A=4sIcpAxo0tofy9s>t| z;f#)h+DlvZ0u&31oe{r9Xwt(3(CkS)?&tG#FRGpgfKvPR?1YTqj1B95)YzZzGR*&% zx%93-zQ+wft8Hi5n{XD@TzY;9aiYK2!9@7pv%RQ^U(vD_sBdksG(@D4MdHCOg!&;L3_+qY*qz}RV3r?j zf-9yPq@b>RkOQ?nw@4mH9&B9%YyDIQN>1SUGSVT?s8O7i-BEk#)bEIMfE^gl0rqEj zBg7QCHv(+`{iG-NgLBj_Z`n(#LkX+U#oOT9;|HJxj!t&MLpkrGaOe@$oCThPkcvPe zC?eud{(BzzvnoQUK$WoU<9>l~I+0kZ#r-@rEMlx0bOcWxj8A-eewi+9@Wg&J%5={$ z7`C}{@n1Q=<5$YL3MkziwD^@;}t*)60@s1Z8`rhwGO{6t5i_#r| z4uK~sbSD9_wTTcs4?f##A}#MPx53xr+plQ4JDmJ_5qJS@(2WM{GmHjSWEc&!2bu=& z02&Mf!`N>RcN3AsI&@s4N_|7^G6MpC8&P}{QKawee)K>x%Ear?9O4fDLiN%%Fen3m zpjkTp0FVrz{gv1OYw&;r`S%B9wxRg*00#&iIt~!-lpuoEYo)khH~N}ij0`}c1p)NK>Y~?l*vg(;q)J;yW?jEMGGI~S=7hm<46ek!6xKQ5dKctRk~JEY zNst}?k?;SgQU=x3T=2EHZF@-vQMHQBPxm6s5Apn%vf$HkO74rc2e2T7kg6#lfRQgjQz!3a+(W7cexpQCpo!%;0s-h+dO^$vb)%9A zok&u@Ql+^280@s&=|xNi0@Hzj&|?4rCZ_`dd6zCB@-E%e|2pJ$36_t)eBon>V3mzn~p$@QufWVfo!ybqx|8H~9LjdxQg#v1&{Y353l?A?-W8LL{P93)+A6n zBV7Ebnh=KR)&o`3t@j(4lU_SNmqw5;XHjemwZnOMgRTG|$xr~ka6gpJ`3|6U2-7$9 zpCb=Az-FK-IxhZ^G95z*H+1!a4t-NWGRUXZ;G=0ODUcK4u3if`6A7 zN^bq$Q~Gb<)4>Rn|Crk?-i`L5!Ny9{6PEQ~(LC=+Tqs|)LyHhjcr@2$engGn zL047=V4m!)@uY)`=W@z&KNACQjJFUy1F^+A3`ux!e1}2bt z`jS6U27&u=i@?H z5FMO@DhMg6#R4by9rIyDpWM=j@}k$#bmo6*C-vrg$O&Od3jSo+SWBi?On);3 zJ!FAn(?izZZ16+FzBJ+h&Tmx)bC(0_9caoQ$vxV0=26mJJ+l$k;-oEYR6Tv3ZBCZ7 z!{$#eTnH0gs_e7;6BqBGFlSp!?Zm1rMB}L&&8n8stP*uHA5W`Z)3J%fTL|T{M>{$h!F@O&T&^~&zs z>On}$4DEuFKkhWFT(CCuN05gicZQ#)(*iXWKMRw5izHm%oo67x{6O3foWV{hJ!Ah# z=9SpJHI8sed4FaJK=5j2@}G1N22E@y4Y_$PCjwWisS3axYLeb9VR!&ANY7{tEQ2I; zBmHm>c=&oeS8JEsi=zOtwV~S{bmq_7b6!BQ5sK+TY9KxEn$&fkj(8SHz-T=_XG)ja zB};U_`Nz?ix2wm!CUub4a^7jIM9LjJql3J8@mpq{8l5|v3+4wbv_k4}LWmv@^lBVL z3INn|!;BsCHcc5cJLlmrI3)%fb=Pa&z-PUhxAAH#lOhW25Tl~E8CFF3?0plw+pZeU z$wDu6ZZsy7UoFQz!)x6cwg=#wYh?FHe%@#U*yMK$!@hCYBWgTe#>*=ESxmF3Jr#B~ zW022ChZ@ptaL!E&I4d#|8_7rg?n2U-}P%U-U5=V1K>PnzASSDLTKSum>_EZ z)-~@>HJn-WBxDV@DM%4r$~xm^Do>5QVb&V(>1|_1q54_EOM^9=qVZw(MAPG%Mbvea zm5q(F3|%?bDDYhcp~jBZ9e|yfv<@V?r_bBET28|*tJdP@Oo9q$k)_y_ChRr@&_pLvb!vm~_mqA;Y2SlEnLVdAL=h2aa8X zIzR~i5Q9_It;;Z$^*+Z$j>2S_>vV+$U2C~m;%az3#F>vtnw5%5aKg?u(9c7@(T8%j zA;G)Tr&ZVZYy&grC%R09<}I@un9qS}OjMIJn_>7xR@*=f6%+;E1P2m5^Veeu#>_YB zgVX~hv*xE=ou<_>wt?Wch()Zog^6VE24{f?2Q5-O2A`UZewfMY&%)Az zq)vn*AHSKjBxx$IP|{smrSf-g~QsT5Z+5EcsJAbpz9QXqX zq3Q9;C3P$hXl`x4*Z0BB34G1%+BJj4Vpezw%-ksp3%rpx2EQ>u{K9sP{ze( z(HW$G1LLC~v}umC?T}sr{ggw(bh1O>me|KDo}q#V^ZQ~ WFDx}mIs;hCQr-oV-~JCfh$n*p literal 0 HcmV?d00001 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/ref_grid.png b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/ref_grid.png new file mode 100644 index 0000000000000000000000000000000000000000..9696f0cd2aeafc9ef053cd481cb4bf68d66a0347 GIT binary patch literal 14414 zcmb7rdpwix|9>Y`-cA)NBr4}a6BuGj1NI$nqSx%F*xfupC7?%TId;D*ID zn|=HC6Zh>qaO3a+;7-U}eIDS4*Td5M+CJ{yUw%tTDsbmWpoL?|zJ2`1_x|_q%X)tj zxXJhEhLsr~JG<}3HDkN5iIp5ku(SD_Yc86V7G6%R<-QS6r`C7< z8Jg0P_iByvl{98^=nSxqW9Z%3!Wg$JzqpLWJXEZ(iGmr5y^hJ_vJ5DKvBEdoR=J!( z)Gk9B#o5wfn;kB=Lggqi5`|4l$KnWwW$a>H_8J3YE;6({)UE(Z#H8!<#BO>8)cq|} z!sPI0cSO;ICnsLmWb42Fb1x_J=KgC7j_qz@2%{xFeH36&%)~A;VgYtoChf)KEVvaK zD?Aq3Q8{w@>zcUCA-I%VQ~Pge^%%n8f=+jYXW_z_H{-mC4V-;&B{l6(96&8$01nUs-A6vryHB3mv^rCpBR6G5`o4+^TQZ0wqSSd!VP= zl`{!4v8xX+OD0En-kts_ul|Gi6MPwvU@~U1@oLP0B$anShX{8ttM{MN!65dwmJxIx zz`Z@S^7ZgQm%)6QgYPCbkD>`Etp^PsbP|X&Zcof}uU^>GiPD!amk3|3dr*nr0(Hj$ZHqomqfEYbiPz3b^Cs+ zjDGc_zUoI6N8`(U)jI$HH~PJLiKtHtag$}yIte(yTt0NK#b0|t@gDkf^xE^(B=PV+ zMrRMfDHmlQasBMA*%6$Dk>gY9CT^n>WMWH(pt^pqw}JPXcmcYqIZnfzD(9eo?F6pl zb~J$1A$Tg|j->f@Ix~{2y@?4KN>IUca9bE0qQk zM^R4Kq(A%54@w*L+;#=}Kb`agCsd^K`ZwwOtr!izdCNIyLU+SfpiLj4jUCMc68-j+ zLR>{-r5|M+Z#Wz6Bp~2RhNAd6vT3oCZ&I?fO_f@@NWX~gi8%Ent)$QF*=j(_Vq&GA z!Tk&@_HweFhJ39Vt`aB0ksUtK%-pK+9E0U-e>G*#joM-6E)M{9)>D^W(>JfAk}#^&M31Ao!3KDv$9L^DP^JoE&CSiw zq_O%j+5!$;*1Zhb#C!b?5NC`=OL1g1%O?zLDfMK{*66gl$*HmW zM@Ohr+V(y4Qw9osHVFWpuc~~X*ItR4N~8x=Xg)fNIYyUuRtW_ZG<>+ZgNh*P6~p|s zh8H>#oR*#e-5-t_o+(ut>{HgT`YN*5+EOZ0;wvEB%`&*JOQ_%rPnA$t$Wqsow;X^6k2~1Q zRPu9HW|5Y{QPt47I&QM<=Hn*CKDL0NRC7j3*87R(jWMoPghA|T9uq(2%Pkn0rkm*j z&+o3p-y3b@LLAn1R`V)Q(qcxYkNVNjmk};q$^d1m{nn{qo9MKY;A)GdEeT~ z%*epwT{A}Fm&!W@OYN%zQm$5sNv|L1Qt>HFiKgStOF^F&ew0_4`z$mr{W5l?ux7xj+VejY7a0NA}+(u)&+X7qmY^5vL&Uww_?53GkyY&dO83p$4; z6#)D_b-^x54o_uPze`N6K^9b;`{J@8wrYeUFNHq6%r3K0wH$a!o5o06*k>wJ)`M7N~bxvy7$GQ@o)T+|oKt6!Fi}AsB7i6*RpJ>og zc~tM>Gbtv^^A|UDuBA?InuB?Ot^3)>+I*Qj92$>-IV19w)xTO{a3u zTH?eUWKFK_^a#N@Hi0&C&%pVlw2q|MBFqTNi@FzIO z2J7C>ZRsT1YnL7}bo9lvdc?KV0@FO>=f1epp0s2FH*k$SIlt!@)DilrHxPl>7huk9 z?pU1YrjDjy^9IB#*|sDsqiFovy`1zCN4ZduQWFfNoh{($HhOZ0C91n~i05HaQD#!Z zY0a+=og2ydIh4`#jfisxI3IR=!2T1>r+pt5cjarRXWcksSAUcg7{)&jHB#Nw`Fn!8 zBZdR;-*~XG!&|SYYlQZw3KYfrR-Q`c9M$OF?J#pL z%%;`1&xE`qJe^GT!TRSJhZ^9R6}*!ka?G7Un`3_(uA(CMLlo;Q}`05oo6aDMyiIJy|n0GXDmeo2t`zT*3 z<>BXUUmVVF^JRWa4_j_rl&0Mtqh6?B<3iz-=3tcHpS!qEanrQjv3k<$s2yp@kljkA zvH!zK-p<|HcCc4*liu~FhqY){bS;$F`WICSSI+*&XF)#m)Nt`*vi^KYKu3Jr*Bfy- z#AECNY?$Dt^Zbp70_+W>!1>@-MSDDd7rJVfk3OH!o_laprz+cJjdw-kWVWc9gvfUm zA*S#A{36n8e2F;0X>(Nf*z+J>xP&MUaq^P{ziStIN`@2xcAvBf=C4^^^WUu$SX)=f(!b3^56BllJ~2H#mioh1Q@RjyP$AWRV}oe@GGSpUQgcK?$9 z`R8V)hwO{`T9g+u%@Yk+`y0rs!vMfKLMg_BF9SCWArBQ~VAnib;VUxoonrbRz$1n}+Iz$lTMUm>1~}6#wA4fx@Mq$*S%g zxOeK;ln38mZ-u5Kql@OE~f@=H8b8)zW zInbB*wo~r&mwZS|dPs&SgLRf#R^5_mSWf=FIGyj%G1m7$@f7$3qESO^Q%9ZZGjZ9} z`73Ge(8oE>OnL2X3lerJ zqq`dXS@V;J$@NDkTr%7^+8?{m%2Q(A%hT)bf@kbwPD5_;YDpfm(V~1{PPHhvJlnG9 zH{(*w=voUfmnwtGsg<=)Zz_1}>(Y?Msk3?encMn&C{4?82}2*&*3ZkRgUzZ(U74XT zR4aIUtAq2>W+q49zPsYUMkI~)K6-RVh4jh4s4Hr{&wNwIy-@o^GoA-?MNomNruTDm z+B_K%SJ63V$Z+ZGuie^-bM#07 z)^Q_EEzR=&kSZ^Zc1`!i3kpByf09Mqf~-~O!^F`;U_Z zAJ+}@I%tk_R~k?Yp&`}v+#&ffe;&N}1|3Y@fm`_UlLH3*W){)1;O^Q0)88CSuXk7V zXo_Bo|8j!$z`{7`0CRx}39gs8wmoEZgdRQm<6zL0uZ=`duyR>v`zNq~HcGYf4WE#6 zUYTmJ{p&2kg?l=cHfRiwE+?7Nxl-9I144TtT3{lPd*E&IsC!QXe8J0djZiJ zwAr!3-SLP}wjM{Tk#|2bopyXQeuV9E!oEju_=|tw2SkvV;|gO{rg?7CJL1lTzox0D z1%L4>4M)&XExEYHQZI_Ct%fB?Z2?w4mP7Qm9ukxJ*teC$ z_#x_+&{Nm_ms2{N0=I}xZP;G+%0c)!YUf5Zy0h>EdgqDB=U&{+8$vc;T!zDYrIc<@W znr`78QSYr&!wVe$p|A?k@cGkO;lAVOi+Q0r9j3pU=P<%1#0*fqlBEwf>{JV;Rz`dpNk0(FzC2oSE}S#g5B^Bp6&dwvZsDZo zQT8jNEU)|AJW~74lNxDmUA?Z6_BVW#b{Z)mw+v}1YxtW!U1(+K95s&(-}e3bxQoKN z?j!qYtdwLO3`f2yfsD;yH?h*xQSQ%N=YO;{wAYUF}v~? z={1UTTmHJ|)wyRNJ3Z$XJ*w}rnrs>z%3z|OWj3)&)Z(pE(ZYXN=4(2tXgO!4ryPits)9Xq1NwP zQfARF($z&S`9Sm(UCqm!@4Pbd-;W}5`!VozUC{C8nU7HqEjT9~b~1D%sh+g266(5I zoH8DEi)pHH+i2IYLAObca}=%iFA+~P<@7@)IO9FEO1}CQ2Dp+v9CiVaw}AlmeD~89bF1|s#i&zud4}RJ!SUSNL_5`1&?Lf9 z)}zQ+pymar+9u~lowcdUO{h|{~=pG{~q;@B!s^BCCeyplj*6D2e>mK^N zvwF<_j{vmy`Zs-~AydNpn0>U1mBYhiY&DE?GjpPOhu3P%?$Ks3?^Pi+j)5ynXqDcn zezlvAp}(e6F(EbKeb`6q6EqbYiz}du2VvC5xs#XwK39m8njF#Soi+m}2ysmw@4OdU zGhBpSh0lAx$&T3?vX7Z({re^^GDh|TVDbV9p#f)3YN_DAjQa1kzjk18+EynpfNLTc zxobM!2`Zllh3gcjed?&g&<@l&y4f*)AoZ5-+7f986x~r`%EgjhSNQ5~N>zE2v!jeqUMye5I0@ zPbR+su?V)#QKGQY%8Pk?IayRT2L=5?x8y8f5yzuvq$R-4XHMj07!413cSdn_e%@a! zccl83s+H`!*LyQ_dlvj2H)E-jE*|3?Rex{RZO}-2B0FcE5U!f%xGjy7RfQ{40y(b2Cn&!{yaOV# zkIDOGd6nWBwiMq1vEdZmZMMJ!s_Gq8%Rw>Qg4}*He3{hTOjIVM z9@PrmN?)$8p54e~P0)IL1*gyj90MVhuoiqv(%cFTBnLqwdX?NK6()VTI;wrGIv~ze zq3*zaa{~zx8141o6k3sEa3n>{m^c<4dkMef5N4AY;UD?j(8@Gg?!~77B zq;L$5dzV(e&GDeeVoktXv;6$ZsB^*cieWLi_RVAY&X!vZR}EY-b?TOle#VXTXqO(> z3;OS)K9#pjv0e~B1?}gmR?j^3qZvd0CC8u+1<}O;`&%gr`a9n4Xy8GAss~Iv1SmRe z)vuMT8m9+nbnf(D$rv*=YlvX}LYqwmnOEoXUWq^C^b2Il*P@D9gaz=gy z?^_elp~m^#Y97Zj3@+?_QPVN{prq|Q`0CP&qT8pDTR*uk%Nenq+lpb*CawYxx3sJ6 zzQ4?Z3n}AyQ9cW+A6+Q&*-=dtd zTWex(%^ZM1Rfsg@NGSIr*=M*%(f_e9QtLa)+a6sIW9b-ScR58B%v><*2IAflA3G8{{(dtJ(dpgM3+Iy(Wo+95tvJd9uPa98!i6r!KPyM{0HB zR`Nh*F8^dTI5M-{_OCl^h*hRJ;+eY@@CJ*hDZ=WHKnvx_K{DRvvca~bbI&xE;2;RC zVA0MGdo_GUYgVPdoZU}byXuZMcHhQ#OV8@E zT!bu++1$-F+W=+zMG7h_3rze%jYW7S^w7_tx>W6F!F{W*2&kkqL2+J9{JqMLaPVP*?v>=|*hocud%bxv%FEj`O?n zEQ&!#{7(jP;mC4lQlNFo1-@DS%jZTO)hHVQEI&{uU`v9N?5Ag>hk`UC+p?PeS8Xxs|~h+TS{R_6qBi8#dyi zlm1ww2xPp&YDyWw&B^!F!X#_bmFhw0uggb98_Qt3=Uh>jDwn)lc>H&0$a-QvkEK-o z`-ZqRors3wZ7Z^Ga3V3Z2%OUF*)$4rIYY+t#k!y@i>o7a7(J_mk9;t#Pv(Vj_pmb4 z90Bl_>c68qRBBThq9NkJVTvhdKq*YLRx!w`KO4yv+sVjwCVFp4K7sL~yermuiMT93 zy9EpKZx%L4fbTOrQdKqXo~C>xOYX0%ej6g4Uk=^FjsA(gf&iP2I*Uw}Q}SWn+*4;c5QaOIIe?G(9E?^(6J{x!DZ1_ls5%4 zgKOt|qu?E1!37-eAdZ2+;PkV)XOAtNU)?Hf_Ygo^@t;xW=J0R#k^f5n3Mzj0(Y^hh z=}2YTmJDoQ7a|RVx`;ol-nfA5C97KU`ST8`YX2MGy64)zq)FQU))w0Om4-|I2psxk z!M{>+f;8RCTE;U?L#ELK9D{4(I_ELBJ0ti!Vg7B#kHMl#Cmv2e_>*rT0acav8;HY1 z?9l!t=Mz|@?nwL1C9(xSW*0ike@>?qyvy!Y~Ha4y_) z99_Okxd-=N%^NTAf9v9aTU4J5?ao6v^`ocmX~*SAvB;7nY0;?dkwC5=S`sgcKe|C- zu18THru}SKFl3#w_6UVv)RKCHn*L459SO0ptRK2vnCzFKj{xJPSwi^s(sRgjvT^#! zV1dh;A1R_uC6aIF7Cgk9&nzrJPDo!n+U-`ld5EDjfGU~Rh>{SGyVn#Xm%1u1I^kFW zF)!11Tv3=rUDmg;FyI*6d4iePq5VA%%8wWDb#*bk|B*Vmq$ino0&P_mcZC`iT3?hY52*Fq%mZ=*DOp=V15fuW zNBE8_%X7DaO2vX-ocK{(lvH`OMt7Q)FNgWP*>I=+o#@9>UiW~TFvg3A5I)XT9~^IZ zh>dCk{^sZDm_%;GXX3yJ7W-{wdU|cQi_$0u3)#ppKkB&pngz*OSrB4g9j6ysIhKlf zGdE$Mooz%8QowAf>Tcx&FL!&$1EpfF2|i0cGzclg!EgAC-8X8c)CR)($)yyN2po3D zITBjbHzV0gy5qh}O!%dg`jjTWn_!<`xDd+T-u`Qkaw>}}qJo6bQg7KhQ_F*t#YYl1 zS6}BzeslzF_pr)T6>NUD7sS>S`!BP+zl>%Q{)T_d3qV`15m+dqD836F>0#SM5hzTNd*_&V|2?xi_`kMp-O8I_^t$ z6yTdAtxPIKy4;+7V<96|k+0{pvXAmF=S|^f{^CB?y+||yaTkjtDbt?+BGF=^B*jlF z79s?oEYA2-(!sU*4V_Qj=-Uv_)6EXW&Iy!87Iv}*v3rnU_%JnZ(9wF(_>TZ3)FCk* zCX6z`r2YBuH*9iwZ8_Jdj-l}d(ROP+#YErbvMpg(W6F0Pw%8t@-Ia|S%J|n6kB=Z&@V0_+eMPGQHZeYBL_ZX!h^4W=qWx^ROH~Zex zjyN-2lvM5Lhh~BK>nRy|`(voYZ09rd*C$>m_=&px$v|}Ze z_w@=?r{MVTP4x5@z~pQ)y5}4@xO*^ny5&B~CdYI!3wtmeG`h3`!<;cDZ;8vBAMJiG zgV`sq(SD2cnBtSa;ok}EB9%i|IQ)n*u&y@+X`?`}J(Ato zNYs1B%}OZ9FY09H;E?KA>sX5J1N-xg_xmXaIOHZ!vANts-l%jJ8v|L;=dWUkh6&Gh zDJH2my0ZHz`4_F~$CXl3R|YayA8h%UJi#4Pr#SIIy1EztVH9uZ{1y8QNB8mMY}e>4 zYG98x9Ck{uP>Y=1Mnxsg)ZXSyw%FX~`WX(A#IX-{7A zP3^t|6+)&mSg~$0vo|p>oLYRk|CO`ew@-X1 z3kqG9((?CrKESUcf@*b1)6ICJ;u8}R8vqD5*l#rVV8&Z$aHgP}L@&I=< zBkRqbtJ8S9g6wG3;j18vOK-&!Q?h=Ey`F{*?qfbLHR;YZxHAwELe*=8iYhDruzHVu zwE(-bQMMXC#vn#UtX`?0x-h*(vXlyUEp&49*#{^o?3Q$#(!^3fv+mnTUw>%pSlYA0 z4c7L#a%=wmlmLilmgFBVM|9)W4hM24Y`JLN(b702LlASMIUDs4(tKl7S+Nq*kN#kx z)4Ln=Yt-k>q8u#m#%zeQ8f9Q=nFaZRIEj+-myd)?Tz2^}f5|(s*+^|L#CB)b>`_^A z>o#hz;p`LINlwJ)Yh}Lgyf#RdHL4AjN6| z*Ybk_`$%pB<5X_4(%;j|_!!+OQX)a1TY2*)rtI8g?t7)VIJ?lBqBW01A1#%fBV2$^ zk&+0j;aOH^uYUKxtat0@F7u;p+zA+llaUU|(D!b?ytdLxfoRO}4_-=2$ms8NrHH>B zLnmBs>h|+m0I3rO8usB`7|WnfW#cztGR=c;v~tCY#~T*io>h`0&Q| z46VYnL22U*0rBWv@eA#~G|@lB7cKcN@bd4cB<&{jb{RZCgB9gPeeSWW!@~WjdR@QP zc9(8ik?p;|&_&q(@9sPQ#-+khzm#s(Hc-=l?wb9bhw!!WO+^Z!qI_LJLgSiqg(j@X zu-P8ht!9FP9%krq1J7DZ9zjG46NJY0aNVy}4XvtZ{^t1e|8^w?nG`~XI2a7~{8~#a z0VrR}oPUJgd$;g4H=dRVi027EKVUQYgcinuZ&FHIN`zAaeAV$(*5*>nyroOkMdIVjXC%>%pb_4c>BY6qtvZa;4sG#hYNef3s!_QUZGV zf7a>bZ%W`>Z5H*|p_p)x;eE8fec`0Ek=_s0MNaD~Vf8`jL8ST|AWc&5P0V|nNnror zgVB9pxwpzMF@5ypTcUIMz`&VRJeMC?vO&*EcO^27#xPSmzobW#CWaMB9z|+oo}vSL z21n(zdCp0F53G6=nLP*F!3S3tf;<|`HfRT+jBG3NB==2}0dLvu%Yy#83Ua8Bx15R| zSE47RE6{yEavEb4{_2&J-b%&=gwlk5c{7Djfva(Ei1Tz{&#J}rDuf)2NV;_e_GJn6GqByvj;_b{mE~#!(uiTpVDj@ zq>IaRLc-IyT#CkEQK7x9tu5Q?bNT zYf*<#&1ODg#Z5&m>*b%M?Y8wZ>NYf(^>vR1*N!5;zh$>5+JxAbk?bQ4BX$-aa~nt@ zx^TO-54D~Urml^xv%{4PLX>@kSy{GDc`uEE+eo-!>UvL2UM-!X`ev%l)_e6ejEwjaR$nkzBx6of5Z%*S6E zFr4u#r0R-{9#MeNO7Yup%4)YAxO95qp-`*}Am1;!8nA5V`zl_%Q!J$=K| zy`yJh&8lHIeBGhPmBHdr7CxtBT^KOnYu~7F?+EUfJ_=L93h_ zzf%5`W?e|h`q@kVW0-*DPQ6r}oV@^5=^uh~3tr$lye_3=S@pOwbKDf>IUmlG4=`k< zIXfdby%LXYdTwS>9QnAi*RVY?q{^t|qK&Qktm+f>jcHiIdFC8_;qSALRwL*uh~aZK zwkpgI-`P=qx$Bq0O=+)vqaKe&3IdVm4g()xzYjY1g#t`iK3VivV2=GQg>w(%=`h!n ztQxI2(^*wl96v)=8P&o@ICo_)a09Ih7j{6|G``zOSY2Rk4d=a0m6zwv<|I zk;a7!j)Pw>e9}FRF9^Ma7X0**@Db8-lTc<|YhlGxPA`MzR?i;^s@0uJ9AGYcKDIxj zo)!oy$oLx2e@<(DbnR6$VIN*0n}FvS#Le}W);IGsZw-nCKOeSZ*R4=34?-UeKq>aL zYVPf6;IIHtJ{$7Qj=pFngppJL#l#Zkc?FgFH3sw^+ScBMc#SuMr06(880-2}q9Bv? zMSyR9r`$+yu<=Giax1&D6G?+0IQRbqn9breNfGK!tpyw zc4^)?)=qxvOuI=aaH#M!KI0BQ^7~0nMqZl~;)VI#Df6@zklb>n`v#LP{zcq#bM*3s zpZ9FEj=(MBY3uZnV-Np9-aUg#(3VVFp3Cp#YYqo>s@bP*t7@wBTVFzDYf7p>b>u;* zZaSUCYov_mrRd2g>NjUFma@C~@umlEY)!f8V1A4-4)8ZVRqvU?+$#(MjYn=j{CJ~G zB4(1e$j6DykG!jkijojyrgUWJL4MuHF8*w#jYFYdY4YL~-V@Xoykl&%RCB$Xjhwtb zTX~TQ7sjqu%_TN%V%fP?E5Is7UfE}s_mPQbfj}+u9-LO_W3&C?3P;T8?Q^^SCsNP<3oqASn8hk{Gtgi?+E_-oIO+{!w{6p0#cY4i?Oqcd?<<7 z+^5BT9DO8GaLGID`g~&U0A<6iKOU!MqzFfO=2_gEN}Pguj$0O1-!oJu zpgn+f3~h}YQ;BE3QIJ1fFE~{olp~3{KN{QWRuZ15H@zo#%Y2O+zS_7y6S7_QX&NwO zMd#3=;dH_dx)C}#5IDScW@hI1FvqxU1fhcK6Iw!O5K-{|xqr9|P~c!`HpVwT0SrNffP zR3h&tb_8DI-({9G81c9;=&HkONYfW_mjR0WT=(&ZIM%=F|C%FWpQ>l|P&NRM#Q02a zMMVmBS3LjCmYGWQ+C(@Hm=%*?4+aZVkx`R!rki^r6>yVB0q_`mBMK63ux9ppmiFudpC{Ak z*7}v7Nr)<`@@@P#QqxJs2X8>i0L;#KZaPI;{qf`>%N$Xh&9*F~nJ(>}SEEyBlSB#$ zM3sw!(hdA0g8+?vuVq9cu=H}OR1gXY#R;#O3C+;bCS3)DJ|6~lN*CDskRnbqumbLF z2m?MAGBoC9)q5OZOBhAFt+lG*WBafnJDA?5lCtqd+%TE|Z2IC295nPu{bZ1^xY%Ki zk6l?%*5duoOw-NWw&gPW1~;pz#N1YjacRdA;xHm}(UDnpaLsIRhW5;Ewl=?EIZ`mV zezMxAlXD=59G^Ra0;<3>_brhZ7i!#oiHV4{l zUp-ZiCggh3fTNrx4oR*wA{(Fr;=5M{xcl#B*m&fx^03>3#;=?lsA)P z$sy>kwRB2dKd1DGx(SJ`su&3*>R}ZdPL17cEK6r9aR{~c)`EG}kgxrhw0xaCoZOK8>&cz&slh(y#?X|dUcXKkRTM&sY)sSIYPE!1sf?cX zq;RxF_*7!D5bAl_>#+493#KUYAT(BzVOfauZ$o(O@!9p1aIFU&&|NN%q>u^wp|R!c z3pTkH%oy=~|GVa7S%iD8Wd8TWDk+VSznI@40DYTrV9iV+hOlcI|16ENBtM3Bno6|R zMWIdNo0vww4guAWPC%ka&Bogru)8BL_ZK$F+ti`3Rpf14QyRQQv}(YGSJuoP{DY$L zpVD^#B>_}&dgrlnFM*hHnE6+FBBP}fn9{V*UgAlV{Xd+*pslnvB6 zASyn7JnM^i!4-~ch7xM%`lp&l)&$nxm~gI40aaASG*lhpoUjjg5m-M0Ahz^6+V$B2 z-#XQLpkSTd@^LZ$-35&z16pIFKA=aZO}bmolU&o|j7k-Lzyh*1=+v7guP?DHd<2@y z^`AYJFfn9iQwbWwYuNk5ubi(6ks?dI6iseUY;Mef0p(xIju*(%z>7N@*pM=0YCKK9 zxOsB#M%D%vRffD~jXR(s1uxxcC``8<6lcf+X=c}U9zi<2TDVGGI;*B4e7c9yfE9gl zCmjKteF#y*xgJD{_;gaTBLs~(LL!TzJyz*__`{)d>ZKdqn7%uEjrvuiU*ALa;Q_@@ zc2f)@1e;_*XBe`7RXA5m7V|&Yh3id9hd^Eq8M3W(6-Q=m=UY}^D@KYm#S-}B+xglr zmu>_%?LqF+N_n5=k)lMvmG#B(n9w&zb$t>I?6O5K*KssI?q%!;@9`}ROpAL%^#a4F zYUzpaL?!BHfN~vW?bU-AKPwXQp&(Kuns60Za>P=3+kHE^?#t;;5F-Hp0s#(uKpU@e zziX7epCPN6O)!iUJgw`$8SDo<1tOuquB4`_YwaB8B=9tskRyR#EhOhOXvL~y+Yh@C zQTdBQ1^+Iz9*|+Mihk?G+T;I93P%&5Z>sMMPWSWxSRO(dVQzffMMS$wPjjn0?HU?I z*0=Yl0YGT~r5xG#dG|?C(ml>z&~aqfQ}xOrgWpn|_E*fQT1sI?4gt4sM*ShAHryAx zzh7jI3*>VXg=7NXf;S&djs(!T#{F8cD9mTe^B(+Z{S8nP17upxks8Bs3p?UGfDmy6 zy#?=YzH0WJfJv6D>GZRyBk`Dvsn-tYtRo(^fw8wpfWbsO;L6m!1l?}~x0Lt;I+Ky+ z2>LU`TET#u#=c}rO|F4*rh%<@Am|MaV+9QpetW>BN>@2TUx3dz9JpsvHD}oLC&R0l zUq{^ltreJZiXH#;%iPk94qi@Ay!zhR;cRUAE?@M#P}6bXmHQQ_;~9N_r~Ck~5dKlZ zF<7cw%;&u;jvp`pp8h?8D&gw|Gyl`(<@LAmfLMn;H|gbz`DY~XKnKlt1fg}d>Ls?Q z05o=p#4%0Cjaj%w`UV(?LMq{SiGNwuo?pu7a0dJtN|QEx4B1%w;dYe#09R)Gzt%5C z_R!ko&bPS#BuM*4G2-qtXHQ35i;3;KMu%FiEb4dv7TtB4|dkx>_8J?5XiWI_g5 zP)H|K9Mk3xiT~gA@#EU({sSL~l!^brIGjcZnFaEo$%Biigkl$f?ob%VqY-cy6VKg* zv+tA^j_Muw3E(}PMS!Hkh7uZs9I;)1RO-?Uo%9BaE9FnW_Xr!9*61n}Kb zDcf`JhldTwD|4cR=6aEs4B`XS48$b4MbSO=s&JBOpNu)Z*z>?*cFsu{j&lp>;L}5a zl=%N|xS*ffz#S9lhK_zUVCAtxi)>GzHi+HW@?R!Ou@L|pfT0W241dRSff%r7=6fGp z`SZmJ2ah7m-zoWZzil|U*XO_2jxr3?$a{3W{zw453)KGOE-n8x$G|5r2`Cp#fZ9cc zD;1n5V)Awk+m@T@umKQ#MUdL=+Z*aR8#vPf1U(bPCo(EwQ#)3e9d+jhUAz_uH_-&1 zoI?P!(FB0?*qW_fe~E&?!GD0ka=R#U;4a|qBeWBXwg+$R5+#@s40#EBA%mKI;aEMl zC*%5C5WsfN0_^?f`90>;nG!yafqZE4+mZvnXarP01@apR7h}$UzVR>dKT-e#W;cM) zYy!+$_r{g|JDCEkrvUWSfFm^;z;V&TGJi!;mjFL6IKcU_NB2+^KuX4L5W0u#|MX+` z<^ML?DEW0jxhre}V>Da10+pMvNjH1a-x}atI08-Q5PpfBi?!a5%50E{wNAzZRfX_U z`~M~IKWayII8RGFx&fx%Fo8uclf#}%L#RC$Og@MLjO{RS4G4^})&eN;Ze^g%FzLEO z0lQus-(zZP&(i^Z4*;Z94s!(l=4CrN{uor2ppbvr1RXe?HOY2z>ds O?}pj!Yn3K$(f<#RUkt+l literal 0 HcmV?d00001 diff --git a/Quasi_static_Finite_strain_Compressible_Elasticity/doc/tooltip b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/tooltip new file mode 100644 index 0000000..a868391 --- /dev/null +++ b/Quasi_static_Finite_strain_Compressible_Elasticity/doc/tooltip @@ -0,0 +1 @@ +Classical quasi-static finite-strain elasticity solving the Cook membrane problem composed of a compressible Neo-Hookean material \ No newline at end of file -- 2.39.5