From: Jean-Paul Pelteret Date: Tue, 14 Feb 2017 08:44:40 +0000 (+0100) Subject: Quasi-static Finite-Strain Quasi-incompressible Visco-elasticity. X-Git-Url: https://gitweb.dealii.org/cgi-bin/gitweb.cgi?a=commitdiff_plain;h=5e2a04dad13b8fc3b51ffabc66c9b71f262437c7;p=code-gallery.git Quasi-static Finite-Strain Quasi-incompressible Visco-elasticity. In this example, we extend step-44 to parallel and implement a rate-dependent constitutive law. We demonstrate the material behaviour by loading a thin strip with a hole in its centre in tension, and track the displacement of several vertices over time. --- diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/CMakeLists.txt b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/CMakeLists.txt new file mode 100644 index 0000000..bd7d78f --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/CMakeLists.txt @@ -0,0 +1,52 @@ +## +# CMake script for the step- tutorial program: +## + +# Set the name of the project and target: +SET(TARGET "viscoelastic_strip_with_hole") + +# Declare all source files the targest consists of: +SET(TARGET_SRC + ${TARGET}.cc +) + +SET(CLEAN_UP_FILES + # a custom list of globs, e.g. *.log *.vtk + *.vtu *.pvtu *.pvd +) + +# Usually, you will not need to modify anything beyond this point... + +CMAKE_MINIMUM_REQUIRED(VERSION 2.8.8) + +FIND_PACKAGE(deal.II 8.5 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() + +# +# Are all dependencies fullfilled? +# +IF(NOT DEAL_II_WITH_CXX11 OR + NOT DEAL_II_WITH_MPI OR + NOT DEAL_II_WITH_METIS OR + NOT DEAL_II_WITH_TRILINOS) + MESSAGE(FATAL_ERROR " +Error! The deal.II library found at ${DEAL_II_PATH} was not configured with + DEAL_II_WITH_CXX11 = ON + DEAL_II_WITH_MPI = ON + DEAL_II_WITH_METIS = ON + DEAL_II_WITH_TRILINOS = ON +One or all of these are OFF in your installation but are required for this code gallery example." + ) +ENDIF() + +DEAL_II_INITIALIZE_CACHED_VARIABLES() +PROJECT(${TARGET}) +DEAL_II_INVOKE_AUTOPILOT() diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/README.md b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/README.md new file mode 100644 index 0000000..50f6ec8 --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/README.md @@ -0,0 +1,134 @@ +## Overview +Many rubber-like materials are not only near-incompressible in nature, but also exhibit a viscoelastic response (within the tested load and time scales). +In this example, we extend the near-incompressible rate-independent constitutive used in step-44 (which implements three-field quasi-static quasi-incompressible finite elasticity) to one that displays rate-dependent behavior. +It may seem that there is a contradiction of terms here, so lets clarify by saying the the problem remains "quasi-static" in the sense that inertial terms remain insignificant, even though the material response itself is rate-dependent. +This implies that, for these fictitious material parameters, it is assumed that the timescale for material relaxation is much longer than that of elastic wave propagation. + +We've also taken the opportunity to extend the code first shown in step-44 to parallel (the primary reason for this contribution), using `Metis` as a grid partitioner and `Trilinos` for linear algebra. +As a motivation as to why one might still choose to use `Metis` (also associated with parallel::shared::Triangulation in step-18, although this triangulation is not used in this instance) over `p4est` (also associated with parallel::distributed::Triangulation) as a grid partitioner, at this point in time it is not possible to use the `hp` finite-element in conjunction with the distributed grid, meaning that this code could not, for example, be readily extended to the application shown in + +* J-P. V. Pelteret, D. Davydov, A. McBride, D. Vu, P. Steinmann, Computational electro-elasticity and magneto-elasticity for quasi-incompressible media immersed in free space. International Journal for Numerical Methods in Engineering, 2016, 108, 1307-1342. DOI: [10.1002/nme.5254](http://doi.org/10.1002/nme.5254) + +The discerning reader will observe that we've chosen to employ `deal.II`'s built in solvers as opposed to using `Trilinos` solvers. +This is because the system matrices `K_Jp` and `K_pJ`, although block diagonal and well conditioned, and for some reason (perhaps pertaining to the negative definite nature of these blocks, or that the entries are very small in magnitude) `Trilinos` solvers are not sufficiently robust to compute inverse matrix-vector multiplication with. +We do stress, however, that to date **no great attempt has been made by the author to overcome this issue** other than by making an entirely different choice of solver. + +### Finite deformation of a thin strip with a hole. + +Various permutations of the problem of an elastomeric strip with a centered cut-out can be found in the literature for solid mechanics, in particular (but not limited to) that pertaining to + +* incompressible elasticity +* elasto-plasticity +* electro-elasticity +* thermo-elasticity. + +Here we implement another permutation (one that is not necessarily benchmarked elsewhere), simply for demonstration purposes. The basic problem configuration is summarized in the following image. + +![Problem geometry](./doc/geometry.png) + +A thin strip of material with a circular hole is (in `3d`) constrained in the `Z` direction and loaded in the direction of its long edge. +In our implementation, this load is applied to the `+Y` surface and may either be displacement control (a `Dirichlet` condition) or a traction load (a `Neumann` boundary condition). +Due to the symmetry of both the geometry and load, the problem can be simplified by modeling only an eighth of the geometry in `3d` or a quarter in `2d`. +By doing so, it it necessary to then implement symmetry conditions on the surfaces coincident with the `X-Z` and `Y-Z` planes (and the `X-Y` plane in `3d`). +The `+X` surface, and that of the hole itself, remain traction-free. + +In three dimensions, the geometry (and a potential partitioning over 8 processors) looks as follows: + +![3d grid with partitioning](./doc/grid_3d-partitioning.png) + +Note that, for this particular formulation, the two-dimensional case corresponds to neither plane-strain nor plane-stress conditions. + +## Requirements +* Version `8.5.0` or greater of `deal.II` +* `C++11` and `MPI` must be enabled +* The following packages must also be enabled: + * `Metis` + * `Trilinos` + +## 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 in serial mode with +``` +make run +``` +and in parallel (in this case, on `4` processors) with +``` +mpirun -np 4 ./viscoelastic_strip_with_hole +``` + +This program can be run in `2d` or `3d`; this choice can be made by making +the appropriate changes in the `main()` function. + + +## 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. DOI: [10.1002/nme.1620371202](http://doi.org/10.1002/nme.1620371202); +* G.A. Holzapfel (2001), Nonlinear Solid Mechanics. A Continuum Approach for Engineering, John Wiley & Sons. ISBN: [978-0-471-82319-3](http://eu.wiley.com/WileyCDA/WileyTitle/productCd-0471823198.html); +* P. Wriggers (2008), Nonlinear finite element methods, Springer. DOI: [10.1007/978-3-540-71001-1](http://doi.org/10.1007/978-3-540-71001-1); +* T.J.R. Hughes (2000), The Finite Element Method: Linear Static and Dynamic Finite Element Analysis, Dover. ISBN: [978-0486411811](http://store.doverpublications.com/0486411818.html) + +The derivation of the finite-element problem, namely the definition and linearization of the residual and their subsequent discretization are quite lengthy and involved. +However, this is already detailed in step-44, some of the aforementioned literature, as well as +* J-P. V. Pelteret, A computational neuromuscular model of the human upper airway with application to the study of obstructive sleep apnoea. PhD Thesis, University of Cape Town, 2013. [http://hdl.handle.net/11427/9519](http://hdl.handle.net/11427/9519) + +and need not be repeated here. +As for the viscoelastic constitutive law (which satisfies the dissipation inequality through the definition of an internal variable that we denote as `Q` in the code), this is derived and presented in detail in +* C. Linder, M. Tkachuk & C. Miehe, A micromechanically motivated diffusion-based transient network model and its incorporation into finite rubber viscoelasticity. Journal of the Mechanics and Physics of Solids, 2011, 59, 2134-2156. DOI: [10.1016/j.jmps.2011.05.005](http://doi.org/10.1016/j.jmps.2011.05.005) + +In particular, the relevant equations from Linder et al.'s work that are implemented in this work are equations 47, 54 and 56. +Note that the time discretization for the rate-dependent internal variable for this single dissipative mechanism is only first-order accurate. + +## Results + +To begin, here is a comparison of the initial grid for the `2d` version of the problem (mirrored about the x-axis and rotated `90` degrees anti-clockwise) + +![Initial grid](./doc/grid_2d-mirrored.png) + +and of the final, displaced grid after the load has been applied and the material is in a (near-completely) relaxed state. + +![Displaced grid](./doc/grid_2d-final_mirrored.png) + +These results, as well as those that follow, were produced using the following material properties: +* The Poisson ratio is `0.4995` +* The elastic shear modulus is `80MPa` +* The viscoelastic shear modulus is `160MPa` +* The viscoelastic relaxation time is `2.5s` + +while the boundary conditions were configured in the following manner: +* The "driving" boundary condition on the short-side (`+Y` face) is + of the `Neumann` variety +* The applied hydrostatic pressure is `-150Pa` (i.e. a tensile load) +* The load ramp time is `1s` + +The following chart highlights the displacement of several vertices and clearly illustrates the viscoelastic nature of the material. + +![Problem geometry](./doc/results-vertex_displacement.png) + +During the initial phase, the load is applied over a period of time much shorter than the material's characteristic relaxation time. +The material therefore exhibits a very stiff response, and relaxes as the load remains constant for the remainder of the simulation. +This deformation that occurs under constant load is commonly known as creep. + +We've been lazy and stopped the simulation slightly prematurely, but it is clear that the material's displacement is moving asymptotically towards a equilibrium solution. +You can also check what the true resting load state is by removing the dissipative mechanism (setting its shear modulus to zero) and rerunning the simulation with the material then exhibiting rate-independent behavior. +Note that in this case, the number of time step over which the load is applied must likely be increased in order to ensure stability of the nonlinear problem. + +### Animations + +Just in case you found the presented chart a little dry and undigestible, below are a couple of animations that demonstrate the viscoelastic nature of the material in a more visually appealing manner. + +**Animation of the evolution of the displacement field** +![Displacement field](./doc/animation-solution_u.gif) + +**Animation of the evolution of the pressure field** +![Pressure field](./doc/animation-solution_p.gif) diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_p.gif b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_p.gif new file mode 100644 index 0000000..51e0025 Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_p.gif differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_u.gif b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_u.gif new file mode 100644 index 0000000..c754ee8 Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/animation-solution_u.gif differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/author b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/author new file mode 100644 index 0000000..6ce7807 --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/author @@ -0,0 +1 @@ +Jean-Paul Pelteret diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/builds-on b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/builds-on new file mode 100644 index 0000000..6f6a518 --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/builds-on @@ -0,0 +1 @@ +step-18 step-44 diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/dependencies b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/dependencies new file mode 100644 index 0000000..2eb2516 --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/dependencies @@ -0,0 +1,4 @@ +DEAL_II_WITH_CXX11 +DEAL_II_WITH_MPI +DEAL_II_WITH_METIS +DEAL_II_WITH_TRILINOS diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/entry-name b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/entry-name new file mode 100644 index 0000000..d492f21 --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/entry-name @@ -0,0 +1 @@ +Quasi-Static Finite-Strain Quasi-incompressible Visco-elasticity diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/geometry.png b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/geometry.png new file mode 100644 index 0000000..cbd8cdd Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/geometry.png differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-final_mirrored.png b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-final_mirrored.png new file mode 100644 index 0000000..e79e16c Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-final_mirrored.png differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-mirrored.png b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-mirrored.png new file mode 100644 index 0000000..711b5c1 Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_2d-mirrored.png differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_3d-partitioning.png b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_3d-partitioning.png new file mode 100644 index 0000000..444ffe2 Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/grid_3d-partitioning.png differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/results-vertex_displacement.png b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/results-vertex_displacement.png new file mode 100644 index 0000000..db27bba Binary files /dev/null and b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/results-vertex_displacement.png differ diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/tooltip b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/tooltip new file mode 100644 index 0000000..374f17f --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/doc/tooltip @@ -0,0 +1 @@ +Quasi-static finite-strain quasi-incompressible rate-dependent elasticity computing the displacement history of a thin viscoelastic strip with a hole. diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/parameters.prm b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/parameters.prm new file mode 100644 index 0000000..7395acf --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/parameters.prm @@ -0,0 +1,102 @@ +# Listing of Parameters +# --------------------- +subsection Boundary conditions + # Driver boundary condition for the problem + set Driver = Neumann + + # Positive stretch applied length-ways to the strip (Driver = Dirichlet) + set Final stretch = 2.0 + + # Hydrostatic pressure applied to top surface of the strip (Driver = Neumann) + # Some reasonablly interesting choices here are -150 and +50 + set Applied pressure = -150.0 + + # Total time over which the stretch/pressure is ramped up + set Load time = 1.0 +end + +subsection Finite element system + # Displacement system polynomial order + set Polynomial degree = 1 + + # Gauss quadrature order + set Quadrature order = 2 +end + + +subsection Geometry + # Total sample length + set Length = 100.0 + + # Total sample width + set Width = 50.0 + + # Total sample thickness + set Thickness = 5.0 + + # Hole diameter + set Hole diameter = 20.0 + + # A geometric factor affecting the discretisation near the hole + set Hole division fraction = 0.5 + + # Number of initial grid subdivisions in the cross-section + set Number of subdivisions in cross-section = 2 + + # Number of initial grid subdivisions through the thickness + set Number of subdivisions thickness = 6 + + # Global refinement level + set Global refinement = 1 + + # Global grid scaling factor + set Grid scale = 1e-3 +end + + +subsection Linear solver + # Linear solver iterations (multiples of the system matrix size) + set Max iteration multiplier = 1 + + # Linear solver residual (scaled by residual norm) + set Residual = 1e-6 + + # Type of solver used to solve the linear system + set Solver type = cg +end + + +subsection Material properties + # Poisson's ratio + set Poisson's ratio = 0.4995 + + # Elastic shear modulus + set Elastic shear modulus = 80.0e6 + + # Viscous shear modulus + set Viscous shear modulus = 160.0e6 + + # Viscous relaxation time + set Viscous relaxation time = 2.5 +end + + +subsection Nonlinear solver + # Number of Newton-Raphson iterations allowed + set Max iterations Newton-Raphson = 10 + + # Displacement error tolerance + set Tolerance displacement = 1.0e-6 + + # Force residual tolerance + set Tolerance force = 1.0e-9 +end + + +subsection Time + # End time + set End time = 25.0 + + # Time step size + set Time step size = 0.1 +end diff --git a/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/viscoelastic_strip_with_hole.cc b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/viscoelastic_strip_with_hole.cc new file mode 100644 index 0000000..73772de --- /dev/null +++ b/Quasi_static_Finite_strain_Quasi_incompressible_ViscoElasticity/viscoelastic_strip_with_hole.cc @@ -0,0 +1,2336 @@ +/* --------------------------------------------------------------------- + * Copyright (C) 2017 by the deal.II authors and + * Jean-Paul Pelteret + * + * 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, 2017 + */ + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +namespace ViscoElasStripHole +{ + using namespace dealii; + namespace LA = TrilinosWrappers; + namespace Parameters + { + struct BoundaryConditions + { + BoundaryConditions(); + + std::string driver; + double stretch; + double pressure; + double load_time; + + const types::boundary_id boundary_id_minus_X; + const types::boundary_id boundary_id_plus_X; + const types::boundary_id boundary_id_minus_Y; + const types::boundary_id boundary_id_plus_Y; + const types::boundary_id boundary_id_minus_Z; + const types::boundary_id boundary_id_plus_Z; + const types::boundary_id boundary_id_hole; + const types::manifold_id manifold_id_hole; + + static void + declare_parameters(ParameterHandler &prm); + void + parse_parameters(ParameterHandler &prm); + }; + BoundaryConditions::BoundaryConditions() + : + driver ("Neumann"), + stretch (2.0), + pressure(0.0), + load_time(2.5), + boundary_id_minus_X (1), + boundary_id_plus_X (2), + boundary_id_minus_Y (3), + boundary_id_plus_Y (4), + boundary_id_minus_Z (5), + boundary_id_plus_Z (6), + boundary_id_hole (10), + manifold_id_hole (10) + { } + void BoundaryConditions::declare_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Boundary conditions"); + { + prm.declare_entry("Driver", "Dirichlet", + Patterns::Selection("Dirichlet|Neumann"), + "Driver boundary condition for the problem"); + prm.declare_entry("Final stretch", "2.0", + Patterns::Double(1.0), + "Positive stretch applied length-ways to the strip"); + prm.declare_entry("Applied pressure", "0.0", + Patterns::Double(-1e3,1e3), + "Hydrostatic pressure applied (in the referential configuration) to the interior surface of the hole"); + prm.declare_entry("Load time", "2.5", + Patterns::Double(0.0), + "Total time over which the stretch/pressure is ramped up"); + } + prm.leave_subsection(); + } + void BoundaryConditions::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Boundary conditions"); + { + driver = prm.get("Driver"); + stretch = prm.get_double("Final stretch"); + pressure = prm.get_double("Applied pressure"); + load_time = prm.get_double("Load time"); + } + prm.leave_subsection(); + } + 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(); + } + struct Geometry + { + double length; + double width; + double thickness; + double hole_diameter; + double hole_division_fraction; + unsigned int n_repetitions_xy; + unsigned int n_repetitions_z; + unsigned int global_refinement; + 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("Length", "100.0", + Patterns::Double(0.0), + "Total sample length"); + prm.declare_entry("Width", "50.0", + Patterns::Double(0.0), + "Total sample width"); + prm.declare_entry("Thickness", "5.0", + Patterns::Double(0.0), + "Total sample thickness"); + prm.declare_entry("Hole diameter", "20.0", + Patterns::Double(0.0), + "Hole diameter"); + prm.declare_entry("Hole division fraction", "0.5", + Patterns::Double(0.0,1.0), + "A geometric factor affecting the discretisation near the hole"); + prm.declare_entry("Number of subdivisions in cross-section", "2", + Patterns::Integer(1.0), + "A factor defining the number of initial grid subdivisions in the cross-section"); + prm.declare_entry("Number of subdivisions thickness", "6", + Patterns::Integer(1.0), + "A factor defining the number of initial grid subdivisions through the thickness"); + prm.declare_entry("Global refinement", "2", + Patterns::Integer(0), + "Global refinement level"); + 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"); + { + length = prm.get_double("Length"); + width = prm.get_double("Width"); + thickness = prm.get_double("Thickness"); + hole_diameter = prm.get_double("Hole diameter"); + hole_division_fraction = prm.get_double("Hole division fraction"); + n_repetitions_xy = prm.get_integer("Number of subdivisions in cross-section"); + n_repetitions_z = prm.get_integer("Number of subdivisions thickness"); + global_refinement = prm.get_integer("Global refinement"); + scale = prm.get_double("Grid scale"); + } + prm.leave_subsection(); + } + struct Materials + { + double nu_e; + double mu_e; + double mu_v; + double tau_v; + 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("Elastic shear modulus", "80.194e6", + Patterns::Double(0.0), + "Elastic shear modulus"); + prm.declare_entry("Viscous shear modulus", "80.194e6", + Patterns::Double(0.0), + "Viscous shear modulus"); + prm.declare_entry("Viscous relaxation time", "2.0", + Patterns::Double(0.0), + "Viscous relaxation time"); + } + prm.leave_subsection(); + } + void Materials::parse_parameters(ParameterHandler &prm) + { + prm.enter_subsection("Material properties"); + { + nu_e = prm.get_double("Poisson's ratio"); + mu_e = prm.get_double("Elastic shear modulus"); + mu_v = prm.get_double("Viscous shear modulus"); + tau_v = prm.get_double("Viscous relaxation time"); + } + prm.leave_subsection(); + } + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + 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(SolverSelector::get_solver_names()), + "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(1.0), + "Linear solver iterations (multiples of the system matrix size)"); + } + 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"); + } + prm.leave_subsection(); + } + 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 displacement", "1.0e-6", + Patterns::Double(0.0), + "Displacement error tolerance"); + prm.declare_entry("Tolerance force", "1.0e-9", + Patterns::Double(0.0), + "Force residual 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(); + } + 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(); + } + struct AllParameters + : public BoundaryConditions, + 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.parse_input(input_file); + parse_parameters(prm); + } + void AllParameters::declare_parameters(ParameterHandler &prm) + { + BoundaryConditions::declare_parameters(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) + { + BoundaryConditions::parse_parameters(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); + } + } + 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; + }; + template + class Material_Compressible_Three_Field_Linear_Viscoelastic + { + public: + Material_Compressible_Three_Field_Linear_Viscoelastic(const double mu_e, + const double nu_e, + const double mu_v, + const double tau_v, + const Time &time) + : + kappa((2.0 * mu_e * (1.0 + nu_e)) / (3.0 * (1.0 - 2.0 * nu_e))), + mu_e(mu_e), + mu_v(mu_v), + tau_v(tau_v), + time(time), + Q_n_t(Physics::Elasticity::StandardTensors::I), + Q_t1(Physics::Elasticity::StandardTensors::I) + { + Assert(kappa > 0, ExcInternalError()); + } + ~Material_Compressible_Three_Field_Linear_Viscoelastic() + {} + + SymmetricTensor<2,dim> + get_tau(const Tensor<2,dim> &F, + const double &p_tilde) const + { + return get_tau_iso(F) + get_tau_vol(F,p_tilde); + } + SymmetricTensor<4,dim> get_Jc(const Tensor<2,dim> &F, + const double &p_tilde) const + { + return get_Jc_iso(F) + get_Jc_vol(F,p_tilde); + } + double + get_dPsi_vol_dJ(const double &J_tilde) const + { + return (kappa / 2.0) * (J_tilde - 1.0 / J_tilde); + } + double + get_d2Psi_vol_dJ2(const double &J_tilde) const + { + return ( (kappa / 2.0) * (1.0 + 1.0 / (J_tilde * J_tilde))); + } + void + update_internal_equilibrium(const Tensor<2, dim> &F, + const double &p_tilde, + const double &J_tilde) + { + const double det_F = determinant(F); + const SymmetricTensor<2,dim> C_bar = std::pow(det_F, -2.0 / dim) * Physics::Elasticity::Kinematics::C(F); + // Linder2011 eq 54 + // Assumes first-oder backward Euler time discretisation + Q_n_t = (1.0/(1.0 + time.get_delta_t()/tau_v))*(Q_t1 + (time.get_delta_t()/tau_v)*invert(C_bar)); + } + void + update_end_timestep() + { + Q_t1 = Q_n_t; + } + + protected: + const double kappa; + const double mu_e; + const double mu_v; + const double tau_v; + const Time &time; + SymmetricTensor<2,dim> Q_n_t; // Value of internal variable at this Newton step and timestep + SymmetricTensor<2,dim> Q_t1; // Value of internal variable at the previous timestep + + SymmetricTensor<2, dim> + get_tau_vol(const Tensor<2,dim> &F, + const double &p_tilde) const + { + const double det_F = determinant(F); + + return p_tilde * det_F * Physics::Elasticity::StandardTensors::I; + } + SymmetricTensor<2, dim> + get_tau_iso(const Tensor<2,dim> &F) const + { + return Physics::Elasticity::StandardTensors::dev_P * get_tau_bar(F); + } + SymmetricTensor<2, dim> + get_tau_bar(const Tensor<2,dim> &F) const + { + const double det_F = determinant(F); + const Tensor<2,dim> F_bar = std::pow(det_F, -1.0 / dim) * F; + const SymmetricTensor<2,dim> b_bar = std::pow(det_F, -2.0 / dim) * symmetrize(F * transpose(F)); + // Elastic Neo-Hookean + Linder2011 eq 47 + return mu_e * b_bar + + mu_v * symmetrize(F_bar*static_cast >(Q_n_t)*transpose(F_bar)); + } + SymmetricTensor<4, dim> get_Jc_vol(const Tensor<2,dim> &F, + const double &p_tilde) const + { + const double det_F = determinant(F); + return p_tilde * det_F + * ( Physics::Elasticity::StandardTensors::IxI + - (2.0 * Physics::Elasticity::StandardTensors::S) ); + } + SymmetricTensor<4, dim> get_Jc_iso(const Tensor<2,dim> &F) const + { + const SymmetricTensor<2, dim> tau_bar = get_tau_bar(F); + const SymmetricTensor<2, dim> tau_iso = get_tau_iso(F); + const SymmetricTensor<4, dim> tau_iso_x_I + = outer_product(tau_iso, + Physics::Elasticity::StandardTensors::I); + const SymmetricTensor<4, dim> I_x_tau_iso + = outer_product(Physics::Elasticity::StandardTensors::I, + tau_iso); + const SymmetricTensor<4, dim> c_bar = get_c_bar(F); + return (2.0 / dim) * trace(tau_bar) + * Physics::Elasticity::StandardTensors::dev_P + - (2.0 / dim) * (tau_iso_x_I + I_x_tau_iso) + + Physics::Elasticity::StandardTensors::dev_P * c_bar + * Physics::Elasticity::StandardTensors::dev_P; + } + SymmetricTensor<4, dim> get_c_bar(const Tensor<2,dim> &F) const + { + // Elastic Neo-Hookean + Linder2011 eq 56 + return -2.0*mu_v*((time.get_delta_t()/tau_v)/(1.0 + time.get_delta_t()/tau_v))*Physics::Elasticity::StandardTensors::S; + } + }; + template + class PointHistory + { + public: + PointHistory() + {} + virtual ~PointHistory() + {} + void + setup_lqp (const Parameters::AllParameters ¶meters, + const Time &time) + { + material.reset(new Material_Compressible_Three_Field_Linear_Viscoelastic( + parameters.mu_e, parameters.nu_e, + parameters.mu_v, parameters.tau_v, + time)); + } + + SymmetricTensor<2, dim> + get_tau(const Tensor<2, dim> &F, + const double &p_tilde) const + { + return material->get_tau(F, p_tilde); + } + SymmetricTensor<4, dim> + get_Jc(const Tensor<2, dim> &F, + const double &p_tilde) const + { + return material->get_Jc(F, p_tilde); + } + double + get_dPsi_vol_dJ(const double &J_tilde) const + { + return material->get_dPsi_vol_dJ(J_tilde); + } + double + get_d2Psi_vol_dJ2(const double &J_tilde) const + { + return material->get_d2Psi_vol_dJ2(J_tilde); + } + void + update_internal_equilibrium(const Tensor<2, dim> &F, + const double &p_tilde, + const double &J_tilde) + { + material->update_internal_equilibrium(F,p_tilde,J_tilde); + } + void + update_end_timestep() + { + material->update_end_timestep(); + } + private: + std_cxx11::shared_ptr< Material_Compressible_Three_Field_Linear_Viscoelastic > material; + }; + template + class Solid + { + public: + Solid(const std::string &input_file); + virtual + ~Solid(); + void + run(); + private: + struct PerTaskData_ASM; + struct ScratchData_ASM; + void + make_grid(); + void + make_2d_quarter_plate_with_hole(Triangulation<2> &tria_2d, + const double half_length, + const double half_width, + const double hole_radius, + const unsigned int n_repetitions_xy = 1, + const double hole_division_fraction = 0.25); + void + setup_system(LA::MPI::BlockVector &solution_delta); + void + determine_component_extractors(); + void + assemble_system(const LA::MPI::BlockVector &solution_delta); + void + assemble_system_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) const; + void + copy_local_to_global_system(const PerTaskData_ASM &data); + void + make_constraints(const int &it_nr); + void + setup_qph(); + void + solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta); + std::pair + solve_linear_system(LA::MPI::BlockVector &newton_update); + LA::MPI::BlockVector + get_solution_total(const LA::MPI::BlockVector &solution_delta) const; + void + update_end_timestep(); + void + output_results(const unsigned int timestep, + const double current_time) const; + void + compute_vertex_positions(std::vector &real_time, + std::vector > > &tracked_vertices, + const LA::MPI::BlockVector &solution_total) const; + + // Parallel communication + MPI_Comm mpi_communicator; + const unsigned int n_mpi_processes; + const unsigned int this_mpi_process; + mutable ConditionalOStream pcout; + + Parameters::AllParameters parameters; + Triangulation triangulation; + Time time; + mutable TimerOutput timer; + CellDataStorage::cell_iterator, + PointHistory > quadrature_point_history; + const unsigned int degree; + const FESystem fe; + DoFHandler dof_handler; + const unsigned int dofs_per_cell; + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Scalar p_fe; + const FEValuesExtractors::Scalar J_fe; + static const unsigned int n_blocks = 3; + static const unsigned int n_components = dim + 2; + static const unsigned int first_u_component = 0; + static const unsigned int p_component = dim; + static const unsigned int J_component = dim + 1; + enum + { + u_block = 0, + p_block = 1, + J_block = 2 + }; + // Block data + std::vector block_component; + + // DoF index data + std::vector all_locally_owned_dofs; + IndexSet locally_owned_dofs; + IndexSet locally_relevant_dofs; + std::vector locally_owned_partitioning; + std::vector locally_relevant_partitioning; + std::vector dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector element_indices_J; + const QGauss qf_cell; + const QGauss qf_face; + const unsigned int n_q_points; + const unsigned int n_q_points_f; + ConstraintMatrix constraints; + LA::BlockSparseMatrix tangent_matrix; + LA::MPI::BlockVector system_rhs; + LA::MPI::BlockVector solution_n; + struct Errors + { + Errors() + : + norm(1.0), u(1.0), p(1.0), J(1.0) + {} + void reset() + { + norm = 1.0; + u = 1.0; + p = 1.0; + J = 1.0; + } + void normalise(const Errors &rhs) + { + if (rhs.norm != 0.0) + norm /= rhs.norm; + if (rhs.u != 0.0) + u /= rhs.u; + if (rhs.p != 0.0) + p /= rhs.p; + if (rhs.J != 0.0) + J /= rhs.J; + } + double norm, u, p, J; + }; + Errors error_residual, error_residual_0, error_residual_norm, error_update, + error_update_0, error_update_norm; + void + get_error_residual(Errors &error_residual); + void + get_error_update(const LA::MPI::BlockVector &newton_update, + Errors &error_update); + std::pair > + get_error_dilation(const LA::MPI::BlockVector &solution_total) const; + void + print_conv_header(); + void + print_conv_footer(const LA::MPI::BlockVector &solution_delta); + }; + template + Solid::Solid(const std::string &input_file) + : + mpi_communicator(MPI_COMM_WORLD), + n_mpi_processes (Utilities::MPI::n_mpi_processes(mpi_communicator)), + this_mpi_process (Utilities::MPI::this_mpi_process(mpi_communicator)), + pcout(std::cout, this_mpi_process == 0), + parameters(input_file), + triangulation(Triangulation::maximum_smoothing), + time(parameters.end_time, parameters.delta_t), + timer(mpi_communicator, + pcout, + TimerOutput::summary, + TimerOutput::wall_times), + degree(parameters.poly_degree), + fe(FE_Q(parameters.poly_degree), dim, // displacement + FE_DGPMonomial(parameters.poly_degree - 1), 1, // pressure + FE_DGPMonomial(parameters.poly_degree - 1), 1), // dilatation + dof_handler(triangulation), + dofs_per_cell (fe.dofs_per_cell), + u_fe(first_u_component), + p_fe(p_component), + J_fe(J_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()) + { + Assert(dim==2 || dim==3, ExcMessage("This problem only works in 2 or 3 space dimensions.")); + determine_component_extractors(); + } + template + Solid::~Solid() + { + dof_handler.clear(); + } + template + void Solid::run() + { + LA::MPI::BlockVector solution_delta; + + make_grid(); + setup_system(solution_delta); + { + ConstraintMatrix constraints; + constraints.close(); + const ComponentSelectFunction + J_mask (J_component, n_components); + VectorTools::project (dof_handler, + constraints, + QGauss(degree+2), + J_mask, + solution_n); + } + output_results(time.get_timestep(), time.current()); + time.increment(); + + // Some points for post-processing + std::vector real_time; + real_time.push_back(0); + std::vector > > tracked_vertices (4); + { + Point p; + p[1] = parameters.length/2.0; + tracked_vertices[0].push_back(p*parameters.scale); + } + { + Point p; + p[1] = parameters.hole_diameter/2.0; + tracked_vertices[1].push_back(p*parameters.scale); + } + { + Point p; + p[0] = parameters.hole_diameter/2.0; + tracked_vertices[2].push_back(p*parameters.scale); + } + { + Point p; + p[0] = parameters.width/2.0; + tracked_vertices[3].push_back(p*parameters.scale); + } + + while (time.current() < time.end()+0.01*time.get_delta_t()) + { + solve_nonlinear_timestep(solution_delta); + solution_n += solution_delta; + solution_delta = 0.0; + output_results(time.get_timestep(), time.current()); + compute_vertex_positions(real_time, + tracked_vertices, + get_solution_total(solution_delta)); + update_end_timestep(); + time.increment(); + } + + pcout << "\n\n*** Spatial position history for tracked vertices ***" << std::endl; + for (unsigned int t=0; t + struct Solid::PerTaskData_ASM + { + FullMatrix cell_matrix; + Vector cell_rhs; + std::vector local_dof_indices; + PerTaskData_ASM(const unsigned int dofs_per_cell) + : + cell_matrix(dofs_per_cell, dofs_per_cell), + cell_rhs(dofs_per_cell), + local_dof_indices(dofs_per_cell) + {} + void reset() + { + cell_matrix = 0.0; + cell_rhs = 0.0; + } + }; + template + struct Solid::ScratchData_ASM + { + const LA::MPI::BlockVector &solution_total; + + // Integration helper + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + + // Quadrature point solution + std::vector > solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_J_total; + + // Shape function values and gradients + std::vector > Nx; + std::vector > > grad_Nx; + std::vector > > symm_grad_Nx; + + ScratchData_ASM(const FiniteElement &fe_cell, + const QGauss &qf_cell, const UpdateFlags uf_cell, + const QGauss & qf_face, const UpdateFlags uf_face, + const LA::MPI::BlockVector &solution_total) + : + solution_total (solution_total), + fe_values_ref(fe_cell, qf_cell, uf_cell), + fe_face_values_ref(fe_cell, qf_face, uf_face), + solution_grads_u_total(qf_cell.size()), + solution_values_p_total(qf_cell.size()), + solution_values_J_total(qf_cell.size()), + 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_ASM(const ScratchData_ASM &rhs) + : + solution_total (rhs.solution_total), + 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()), + solution_grads_u_total(rhs.solution_grads_u_total), + solution_values_p_total(rhs.solution_values_p_total), + solution_values_J_total(rhs.solution_values_J_total), + Nx(rhs.Nx), + grad_Nx(rhs.grad_Nx), + symm_grad_Nx(rhs.symm_grad_Nx) + {} + void reset() + { + const unsigned int n_q_points = solution_grads_u_total.size(); + const unsigned int n_dofs_per_cell = Nx[0].size(); + + Assert(solution_grads_u_total.size() == n_q_points, + ExcInternalError()); + Assert(solution_values_p_total.size() == n_q_points, + ExcInternalError()); + Assert(solution_values_J_total.size() == n_q_points, + ExcInternalError()); + Assert(Nx.size() == n_q_points, + ExcInternalError()); + Assert(grad_Nx.size() == n_q_points, + ExcInternalError()); + Assert(symm_grad_Nx.size() == n_q_points, + ExcInternalError()); + + 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()); + + solution_grads_u_total[q_point] = 0.0; + solution_values_p_total[q_point] = 0.0; + solution_values_J_total[q_point] = 0.0; + 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; + } + } + } + }; + template <> + void Solid<2>::make_grid() + { + const int dim = 2; + const double tol = 1e-12; + make_2d_quarter_plate_with_hole(triangulation, + parameters.length/2.0, + parameters.width/2.0, + parameters.hole_diameter/2.0, + parameters.n_repetitions_xy, + parameters.hole_division_fraction); + + // Clear boundary ID's + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + cell->face(face)->set_all_boundary_ids(0); + } + } + + // Set boundary IDs and and manifolds + const Point centre (0,0); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + // Set boundary IDs + if (std::abs(cell->face(face)->center()[0] - 0.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X); + } + else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X); + } + else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y); + } + else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y); + } + else + { + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_hole); + break; + } + } + + // Set manifold IDs + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + if (std::abs(cell->vertex(vertex).distance(centre) - parameters.hole_diameter/2.0) < tol) + { + cell->face(face)->set_manifold_id(parameters.manifold_id_hole); + break; + } + } + } + static SphericalManifold spherical_manifold (centre); + triangulation.set_manifold(parameters.manifold_id_hole,spherical_manifold); + triangulation.refine_global(parameters.global_refinement); + GridTools::scale(parameters.scale,triangulation); + } + template <> + void Solid<3>::make_grid() + { + const int dim = 3; + const double tol = 1e-12; + Triangulation<2> tria_2d; + make_2d_quarter_plate_with_hole(tria_2d, + parameters.length/2.0, + parameters.width/2.0, + parameters.hole_diameter/2.0, + parameters.n_repetitions_xy, + parameters.hole_division_fraction); + GridGenerator::extrude_triangulation(tria_2d, + parameters.n_repetitions_z+1, + parameters.thickness/2.0, + triangulation); + + // Clear boundary ID's + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + cell->face(face)->set_all_boundary_ids(0); + } + } + + // Set boundary IDs and and manifolds + const Point direction (0,0,1); + const Point centre (0,0,0); + for (typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(); + cell != triangulation.end(); ++cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + { + // Set boundary IDs + if (std::abs(cell->face(face)->center()[0] - 0.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_minus_X); + } + else if (std::abs(cell->face(face)->center()[0] - parameters.width/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_plus_X); + } + else if (std::abs(cell->face(face)->center()[1] - 0.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Y); + } + else if (std::abs(cell->face(face)->center()[1] - parameters.length/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Y); + } + else if (std::abs(cell->face(face)->center()[2] - 0.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_minus_Z); + } + else if (std::abs(cell->face(face)->center()[2] - parameters.thickness/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_plus_Z); + } + else + { + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + { + // Project the cell vertex to the XY plane and + // test the distance from the cylinder axis + Point vertex_proj = cell->vertex(vertex); + vertex_proj[2] = 0.0; + if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < tol) + { + cell->face(face)->set_boundary_id(parameters.boundary_id_hole); + break; + } + } + } + + // Set manifold IDs + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + { + // Project the cell vertex to the XY plane and + // test the distance from the cylinder axis + Point vertex_proj = cell->vertex(vertex); + vertex_proj[2] = 0.0; + if (std::abs(vertex_proj.distance(centre) - parameters.hole_diameter/2.0) < 1e-12) + { + // Set manifold ID on face and edges + cell->face(face)->set_all_manifold_ids(parameters.manifold_id_hole); + break; + } + } + } + } + static CylindricalManifold cylindrical_manifold (direction,centre); + triangulation.set_manifold(parameters.manifold_id_hole,cylindrical_manifold); + triangulation.refine_global(parameters.global_refinement); + GridTools::scale(parameters.scale,triangulation); + } + template + void Solid::make_2d_quarter_plate_with_hole(Triangulation<2> &tria_2d, + const double half_length, + const double half_width, + const double hole_radius, + const unsigned int n_repetitions_xy, + const double hole_division_fraction) + { + const double length = 2.0*half_length; + const double width = 2.0*half_width; + const double hole_diameter = 2.0*hole_radius; + + const double internal_width = hole_diameter + hole_division_fraction*(width - hole_diameter); + Triangulation<2> tria_quarter_plate_hole; + { + Triangulation<2> tria_plate_hole; + GridGenerator::hyper_cube_with_cylindrical_hole (tria_plate_hole, + hole_diameter/2.0, + internal_width/2.0); + + std::set::active_cell_iterator > cells_to_remove; + for (typename Triangulation<2>::active_cell_iterator + cell = tria_plate_hole.begin_active(); + cell != tria_plate_hole.end(); ++cell) + { + // Remove all cells that are not in the first quadrant + if (cell->center()[0] < 0.0 || cell->center()[1] < 0.0) + cells_to_remove.insert(cell); + } + Assert(cells_to_remove.size() > 0, ExcInternalError()); + Assert(cells_to_remove.size() != tria_plate_hole.n_active_cells(), ExcInternalError()); + GridGenerator::create_triangulation_with_removed_cells(tria_plate_hole,cells_to_remove,tria_quarter_plate_hole); + } + + Triangulation<2> tria_cut_plate; + { + Triangulation<2> tria_plate; + // Subdivide the plate so that we're left one + // cell to remove (we'll replace this with the + // plate with the hole) and then make the + // rest of the subdivisions so that we're left + // with cells with a decent aspect ratio + std::vector > step_sizes; + { + std::vector subdivision_width; + subdivision_width.push_back(internal_width/2.0); + const double width_remaining = (width - internal_width)/2.0; + const unsigned int n_subs = std::max(1.0,std::ceil(width_remaining/(internal_width/2.0))); + Assert(n_subs>0, ExcInternalError()); + for (unsigned int s=0; s subdivision_length; + subdivision_length.push_back(internal_width/2.0); + const double length_remaining = (length - internal_width)/2.0; + const unsigned int n_subs = std::max(1.0,std::ceil(length_remaining/(internal_width/2.0))); + Assert(n_subs>0, ExcInternalError()); + for (unsigned int s=0; s(0.0, 0.0), + Point<2>(width/2.0, length/2.0)); + + std::set::active_cell_iterator > cells_to_remove; + for (typename Triangulation<2>::active_cell_iterator + cell = tria_plate.begin_active(); + cell != tria_plate.end(); ++cell) + { + // Remove all cells that are in the first quadrant + if (cell->center()[0] < internal_width/2.0 && cell->center()[1] < internal_width/2.0) + cells_to_remove.insert(cell); + } + Assert(cells_to_remove.size() > 0, ExcInternalError()); + Assert(cells_to_remove.size() != tria_plate.n_active_cells(), ExcInternalError()); + GridGenerator::create_triangulation_with_removed_cells(tria_plate,cells_to_remove,tria_cut_plate); + } + + Triangulation<2> tria_2d_not_flat; + GridGenerator::merge_triangulations(tria_quarter_plate_hole, + tria_cut_plate, + tria_2d_not_flat); + + // Attach a manifold to the curved boundary and refine + // Note: We can only guarentee that the vertices sit on + // the curve, so we must test with their position instead + // of the cell centre. + const Point<2> centre_2d (0,0); + for (typename Triangulation<2>::active_cell_iterator + cell = tria_2d_not_flat.begin_active(); + cell != tria_2d_not_flat.end(); ++cell) + { + for (unsigned int face=0; face::faces_per_cell; ++face) + if (cell->face(face)->at_boundary()) + for (unsigned int vertex=0; vertex::vertices_per_face; ++vertex) + if (std::abs(cell->vertex(vertex).distance(centre_2d) - hole_diameter/2.0) < 1e-12) + { + cell->face(face)->set_manifold_id(10); + break; + } + } + SphericalManifold<2> spherical_manifold_2d (centre_2d); + tria_2d_not_flat.set_manifold(10,spherical_manifold_2d); + tria_2d_not_flat.refine_global(std::max (1U, n_repetitions_xy)); + tria_2d_not_flat.set_manifold(10); // Clear manifold + + GridGenerator::flatten_triangulation(tria_2d_not_flat,tria_2d); + } + template + void Solid::setup_system(LA::MPI::BlockVector &solution_delta) + { + timer.enter_subsection("Setup system"); + pcout << "Setting up linear system..." << std::endl; + + // Partition triangulation + GridTools::partition_triangulation (n_mpi_processes, + triangulation); + + block_component = std::vector (n_components, u_block); // Displacement + block_component[p_component] = p_block; // Pressure + block_component[J_component] = J_block; // Dilatation + dof_handler.distribute_dofs(fe); + DoFRenumbering::Cuthill_McKee(dof_handler); + DoFRenumbering::component_wise(dof_handler, block_component); + + // Count DoFs in each block + dofs_per_block.clear(); + dofs_per_block.resize(n_blocks); + DoFTools::count_dofs_per_block(dof_handler, dofs_per_block, + block_component); + + all_locally_owned_dofs = DoFTools::locally_owned_dofs_per_subdomain (dof_handler); + std::vector all_locally_relevant_dofs + = DoFTools::locally_relevant_dofs_per_subdomain (dof_handler); + + locally_owned_dofs.clear(); + locally_owned_partitioning.clear(); + Assert(all_locally_owned_dofs.size() > this_mpi_process, ExcInternalError()); + locally_owned_dofs = all_locally_owned_dofs[this_mpi_process]; + + locally_relevant_dofs.clear(); + locally_relevant_partitioning.clear(); + Assert(all_locally_relevant_dofs.size() > this_mpi_process, ExcInternalError()); + locally_relevant_dofs = all_locally_relevant_dofs[this_mpi_process]; + + locally_owned_partitioning.reserve(n_blocks); + locally_relevant_partitioning.reserve(n_blocks); + for (unsigned int b=0; b coupling(n_components, n_components); + for (unsigned int ii = 0; ii < n_components; ++ii) + for (unsigned int jj = 0; jj < n_components; ++jj) + if (((ii < p_component) && (jj == J_component)) + || ((ii == J_component) && (jj < p_component)) + || ((ii == p_component) && (jj == p_component))) + coupling[ii][jj] = DoFTools::none; + else + coupling[ii][jj] = DoFTools::always; + + TrilinosWrappers::BlockSparsityPattern bsp (locally_owned_partitioning, + locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator); + DoFTools::make_sparsity_pattern (dof_handler, bsp, + constraints, false, + this_mpi_process); + bsp.compress(); + tangent_matrix.reinit (bsp); + + // We then set up storage vectors + system_rhs.reinit(locally_owned_partitioning, + mpi_communicator); + solution_n.reinit(locally_owned_partitioning, + mpi_communicator); + solution_delta.reinit(locally_owned_partitioning, + mpi_communicator); + setup_qph(); + timer.leave_subsection(); + } + template + void + Solid::determine_component_extractors() + { + element_indices_u.clear(); + element_indices_p.clear(); + element_indices_J.clear(); + for (unsigned int k = 0; k < fe.dofs_per_cell; ++k) + { + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_block) + element_indices_u.push_back(k); + else if (k_group == p_block) + element_indices_p.push_back(k); + else if (k_group == J_block) + element_indices_J.push_back(k); + else + { + Assert(k_group <= J_block, ExcInternalError()); + } + } + } + template + void Solid::setup_qph() + { + pcout << "Setting up quadrature point data..." << std::endl; + quadrature_point_history.initialize(triangulation.begin_active(), + triangulation.end(), + n_q_points); + FilteredIterator::active_cell_iterator> + cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.begin_active()), + endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.end()); + for (; cell!=endc; ++cell) + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->setup_lqp(parameters, time); + } + } + template + void + Solid::solve_nonlinear_timestep(LA::MPI::BlockVector &solution_delta) + { + pcout << std::endl + << "Timestep " << time.get_timestep() << " @ " + << time.current() << "s of " + << time.end() << "s" << std::endl; + LA::MPI::BlockVector newton_update(locally_owned_partitioning, + mpi_communicator); + 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(); + unsigned int newton_iteration = 0; + for (; newton_iteration < parameters.max_iterations_NR; + ++newton_iteration) + { + pcout << " " << std::setw(2) << newton_iteration << " " << std::flush; + make_constraints(newton_iteration); + assemble_system(solution_delta); + get_error_residual(error_residual); + if (newton_iteration == 0) + error_residual_0 = error_residual; + 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) ) + { + pcout << " CONVERGED! " << std::endl; + print_conv_footer(solution_delta); + break; + } + 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; + error_update_norm = error_update; + error_update_norm.normalise(error_update_0); + solution_delta += newton_update; + newton_update = 0.0; + pcout << " | " << 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_residual_norm.p << " " << error_residual_norm.J + << " " << error_update_norm.norm << " " << error_update_norm.u + << " " << error_update_norm.p << " " << error_update_norm.J + << " " << std::endl; + } + AssertThrow (newton_iteration <= parameters.max_iterations_NR, + ExcMessage("No convergence in nonlinear solver!")); + } + template + void Solid::print_conv_header() + { + pcout << std::string(132,'_') << std::endl; + pcout << " SOLVER STEP " + << " | LIN_IT LIN_RES RES_NORM " + << " RES_U RES_P RES_J NU_NORM " + << " NU_U NU_P NU_J " << std::endl; + pcout << std::string(132,'_') << std::endl; + } + template + void Solid::print_conv_footer(const LA::MPI::BlockVector &solution_delta) + { + pcout << std::string(132,'_') << std::endl; + const std::pair > error_dil = get_error_dilation(get_solution_total(solution_delta)); + pcout << "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 + << "Dilatation:\t" << error_dil.first << std::endl + << "v / V_0:\t" << error_dil.second.second << " / " << error_dil.second.first + << " = " << (error_dil.second.second/error_dil.second.first) << std::endl; + } + template + std::pair > + Solid::get_error_dilation(const LA::MPI::BlockVector &solution_total) const + { + double vol_reference = 0.0; + double vol_current = 0.0; + double dil_L2_error = 0.0; + FEValues fe_values_ref(fe, qf_cell, + update_values | update_gradients | update_JxW_values); + std::vector > solution_grads_u_total (qf_cell.size()); + std::vector solution_values_J_total (qf_cell.size()); + FilteredIterator::active_cell_iterator> + cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.begin_active()), + endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.end()); + for (; cell != endc; ++cell) + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + fe_values_ref.reinit(cell); + fe_values_ref[u_fe].get_function_gradients(solution_total, + solution_grads_u_total); + fe_values_ref[J_fe].get_function_values(solution_total, + solution_values_J_total); + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const double det_F_qp = determinant(Physics::Elasticity::Kinematics::F(solution_grads_u_total[q_point])); + const double J_tilde_qp = solution_values_J_total[q_point]; + const double the_error_qp_squared = std::pow((det_F_qp - J_tilde_qp), + 2); + const double JxW = fe_values_ref.JxW(q_point); + dil_L2_error += the_error_qp_squared * JxW; + vol_reference += JxW; + vol_current += det_F_qp * JxW; + } + } + Assert(vol_current > 0.0, ExcInternalError()); + // Sum across all porcessors + dil_L2_error = Utilities::MPI::sum(dil_L2_error,mpi_communicator); + vol_reference = Utilities::MPI::sum(vol_reference,mpi_communicator); + vol_current = Utilities::MPI::sum(vol_current,mpi_communicator); + + return std::make_pair(std::sqrt(dil_L2_error), + std::make_pair(vol_reference,vol_current)); + } + template + void Solid::get_error_residual(Errors &error_residual) + { + // Construct a residual vector that has the values for all of its + // constrained DoFs set to zero. + LA::MPI::BlockVector error_res (system_rhs); + constraints.set_zero(error_res); + error_residual.norm = error_res.l2_norm(); + error_residual.u = error_res.block(u_block).l2_norm(); + error_residual.p = error_res.block(p_block).l2_norm(); + error_residual.J = error_res.block(J_block).l2_norm(); + } + template + void Solid::get_error_update(const LA::MPI::BlockVector &newton_update, + Errors &error_update) + { + // Construct a update vector that has the values for all of its + // constrained DoFs set to zero. + LA::MPI::BlockVector error_ud (newton_update); + constraints.set_zero(error_ud); + error_update.norm = error_ud.l2_norm(); + error_update.u = error_ud.block(u_block).l2_norm(); + error_update.p = error_ud.block(p_block).l2_norm(); + error_update.J = error_ud.block(J_block).l2_norm(); + } + template + LA::MPI::BlockVector + Solid::get_solution_total(const LA::MPI::BlockVector &solution_delta) const + { + // Cell interpolation -> Ghosted vector + LA::MPI::BlockVector solution_total (locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator, + /*vector_writable = */ false); + LA::MPI::BlockVector tmp (solution_total); + solution_total = solution_n; + tmp = solution_delta; + solution_total += tmp; + return solution_total; + } + template + void Solid::assemble_system(const LA::MPI::BlockVector &solution_delta) + { + timer.enter_subsection("Assemble system"); + pcout << " ASM_SYS " << std::flush; + tangent_matrix = 0.0; + system_rhs = 0.0; + const LA::MPI::BlockVector solution_total(get_solution_total(solution_delta)); + const UpdateFlags uf_cell(update_values | + update_gradients | + update_JxW_values); + const UpdateFlags uf_face(update_values | + update_normal_vectors | + update_JxW_values); + PerTaskData_ASM per_task_data(dofs_per_cell); + ScratchData_ASM scratch_data(fe, qf_cell, uf_cell, qf_face, uf_face, solution_total); + + FilteredIterator::active_cell_iterator> + cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.begin_active()), + endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.end()); + for (; cell != endc; ++cell) + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + assemble_system_one_cell(cell, scratch_data, per_task_data); + copy_local_to_global_system(per_task_data); + } + tangent_matrix.compress(VectorOperation::add); + system_rhs.compress(VectorOperation::add); + timer.leave_subsection(); + } + template + void Solid::copy_local_to_global_system(const PerTaskData_ASM &data) + { + constraints.distribute_local_to_global(data.cell_matrix, data.cell_rhs, + data.local_dof_indices, + tangent_matrix, system_rhs); + } + template + void + Solid::assemble_system_one_cell(const typename DoFHandler::active_cell_iterator &cell, + ScratchData_ASM &scratch, + PerTaskData_ASM &data) const + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + + data.reset(); + scratch.reset(); + scratch.fe_values_ref.reinit(cell); + cell->get_dof_indices(data.local_dof_indices); + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + + // Update quadrature point solution + scratch.fe_values_ref[u_fe].get_function_gradients(scratch.solution_total, + scratch.solution_grads_u_total); + scratch.fe_values_ref[p_fe].get_function_values(scratch.solution_total, + scratch.solution_values_p_total); + scratch.fe_values_ref[J_fe].get_function_values(scratch.solution_total, + scratch.solution_values_J_total); + + // Update shape functions and their gradients (push-forward) + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]); + const Tensor<2, dim> F_inv = invert(F); + + 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_block) + { + 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 if (k_group == p_block) + scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, + q_point); + else if (k_group == J_block) + scratch.Nx[q_point][k] = scratch.fe_values_ref[J_fe].value(k, + q_point); + else + Assert(k_group <= J_block, ExcInternalError()); + } + } + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + { + const SymmetricTensor<2, dim> &I = Physics::Elasticity::StandardTensors::I; + const Tensor<2, dim> F = Physics::Elasticity::Kinematics::F(scratch.solution_grads_u_total[q_point]); + const double det_F = determinant(F); + const double &p_tilde = scratch.solution_values_p_total[q_point]; + const double &J_tilde = scratch.solution_values_J_total[q_point]; + Assert(det_F > 0, ExcInternalError()); + + { + PointHistory *lqph_q_point_nc = const_cast*>(lqph[q_point].get()); + lqph_q_point_nc->update_internal_equilibrium(F,p_tilde,J_tilde); + } + + const SymmetricTensor<2, dim> tau = lqph[q_point]->get_tau(F,p_tilde); + const Tensor<2, dim> tau_ns (tau); + const SymmetricTensor<4, dim> Jc = lqph[q_point]->get_Jc(F,p_tilde); + const double dPsi_vol_dJ = lqph[q_point]->get_dPsi_vol_dJ(J_tilde); + const double d2Psi_vol_dJ2 = lqph[q_point]->get_d2Psi_vol_dJ2(J_tilde); + + const std::vector &Nx = scratch.Nx[q_point]; + const std::vector > &grad_Nx = scratch.grad_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); + + 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; + if (i_group == u_block) + data.cell_rhs(i) -= (symm_grad_Nx[i] * tau) * JxW; + else if (i_group == p_block) + data.cell_rhs(i) -= Nx[i] * (det_F - J_tilde) * JxW; + else if (i_group == J_block) + data.cell_rhs(i) -= Nx[i] * (dPsi_vol_dJ - p_tilde) * JxW; + else + Assert(i_group <= J_block, ExcInternalError()); + + 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; + if ((i_group == u_block) && (j_group == u_block)) + { + 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_ns + * grad_Nx[j][component_j] * JxW; + } + else if ((i_group == u_block) && (j_group == p_block)) + { + data.cell_matrix(i, j) += (symm_grad_Nx[i] * I) + * Nx[j] * det_F + * JxW; + } + else if ((i_group == p_block) && (j_group == u_block)) + { + data.cell_matrix(i, j) += Nx[i] * det_F + * (symm_grad_Nx[j] * I) + * JxW; + } + else if ((i_group == p_block) && (j_group == J_block)) + data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW; + else if ((i_group == J_block) && (j_group == p_block)) + data.cell_matrix(i, j) -= Nx[i] * Nx[j] * JxW; + else if ((i_group == J_block) && (j_group == J_block)) + data.cell_matrix(i, j) += Nx[i] * d2Psi_vol_dJ2 * Nx[j] * JxW; + else + Assert((i_group <= J_block) && (j_group <= J_block), + ExcInternalError()); + } + } + } + + 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); + + if (parameters.driver == "Neumann") + for (unsigned int face = 0; face < GeometryInfo::faces_per_cell; + ++face) + if (cell->face(face)->at_boundary() == true + && cell->face(face)->boundary_id() == parameters.boundary_id_plus_Y) + { + 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) + { + const Tensor<1, dim> &N = + scratch.fe_face_values_ref.normal_vector(f_q_point); + static const double pressure_nom = parameters.pressure + / (parameters.scale * parameters.scale); + const double time_ramp = (time.current() < parameters.load_time ? + time.current() / parameters.load_time : 1.0); + const double pressure = -pressure_nom * time_ramp; + const Tensor<1, dim> traction = pressure * N; + 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_block) + { + 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; + } + } + } + } + } + template + void Solid::make_constraints(const int &it_nr) + { + pcout << " CST " << std::flush; + if (it_nr > 1) + return; + constraints.clear(); + const bool apply_dirichlet_bc = (it_nr == 0); + const FEValuesExtractors::Scalar x_displacement(0); + const FEValuesExtractors::Scalar y_displacement(1); + { + const int boundary_id = parameters.boundary_id_minus_X; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(x_displacement)); + } + { + const int boundary_id = parameters.boundary_id_minus_Y; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + } + if (dim==3) + { + const FEValuesExtractors::Scalar z_displacement(2); + { + const int boundary_id = parameters.boundary_id_minus_Z; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + { + const int boundary_id = parameters.boundary_id_plus_Z; + if (apply_dirichlet_bc == true) + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(z_displacement)); + } + } + if (parameters.driver == "Dirichlet") + { + const int boundary_id = parameters.boundary_id_plus_Y; + if (apply_dirichlet_bc == true) + { + + if (time.current() < parameters.load_time+0.01*time.get_delta_t()) + { + const double delta_length = parameters.length*(parameters.stretch - 1.0)*parameters.scale; + const unsigned int n_stretch_steps = parameters.load_time/time.get_delta_t(); + const double delta_u_y = delta_length/2.0/n_stretch_steps; + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ConstantFunction(delta_u_y,n_components), + constraints, + fe.component_mask(y_displacement)); + } + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + } + else + VectorTools::interpolate_boundary_values(dof_handler, + boundary_id, + ZeroFunction(n_components), + constraints, + fe.component_mask(y_displacement)); + } + constraints.close(); + } + template + std::pair + Solid::solve_linear_system(LA::MPI::BlockVector &newton_update) + { + unsigned int lin_it = 0; + double lin_res = 0.0; + + timer.enter_subsection("Linear solver"); + pcout << " SLV " << std::flush; + + const LA::MPI::Vector &f_u = system_rhs.block(u_block); + const LA::MPI::Vector &f_p = system_rhs.block(p_block); + const LA::MPI::Vector &f_J = system_rhs.block(J_block); + LA::MPI::Vector &d_u = newton_update.block(u_block); + LA::MPI::Vector &d_p = newton_update.block(p_block); + LA::MPI::Vector &d_J = newton_update.block(J_block); + const auto K_uu = linear_operator(tangent_matrix.block(u_block, u_block)); + const auto K_up = linear_operator(tangent_matrix.block(u_block, p_block)); + const auto K_pu = linear_operator(tangent_matrix.block(p_block, u_block)); + const auto K_Jp = linear_operator(tangent_matrix.block(J_block, p_block)); + const auto K_JJ = linear_operator(tangent_matrix.block(J_block, J_block)); + + LA::PreconditionJacobi preconditioner_K_Jp_inv; + preconditioner_K_Jp_inv.initialize( + tangent_matrix.block(J_block, p_block), + LA::PreconditionJacobi::AdditionalData()); + ReductionControl solver_control_K_Jp_inv ( + tangent_matrix.block(J_block, p_block).m() * parameters.max_iterations_lin, + 1.0e-30, 1e-6); + dealii::SolverCG solver_K_Jp_inv (solver_control_K_Jp_inv); + + const auto K_Jp_inv = inverse_operator(K_Jp, + solver_K_Jp_inv, + preconditioner_K_Jp_inv); + const auto K_pJ_inv = transpose_operator(K_Jp_inv); + const auto K_pp_bar = K_Jp_inv * K_JJ * K_pJ_inv; + const auto K_uu_bar_bar = K_up * K_pp_bar * K_pu; + const auto K_uu_con = K_uu + K_uu_bar_bar; + + LA::PreconditionAMG preconditioner_K_con_inv; + preconditioner_K_con_inv.initialize( + tangent_matrix.block(u_block, u_block), + LA::PreconditionAMG::AdditionalData( + true /*elliptic*/, + (parameters.poly_degree > 1 /*higher_order_elements*/)) ); + ReductionControl solver_control_K_con_inv ( + tangent_matrix.block(u_block, u_block).m() * parameters.max_iterations_lin, + 1.0e-30, parameters.tol_lin); + dealii::SolverSelector solver_K_con_inv; + solver_K_con_inv.select(parameters.type_lin); + solver_K_con_inv.set_control(solver_control_K_con_inv); + const auto K_uu_con_inv = inverse_operator(K_uu_con, + solver_K_con_inv, + preconditioner_K_con_inv); + + d_u = K_uu_con_inv*(f_u - K_up*(K_Jp_inv*f_J - K_pp_bar*f_p)); + lin_it = solver_control_K_con_inv.last_step(); + lin_res = solver_control_K_con_inv.last_value(); + timer.leave_subsection(); + + timer.enter_subsection("Linear solver postprocessing"); + d_J = K_pJ_inv*(f_p - K_pu*d_u); + d_p = K_Jp_inv*(f_J - K_JJ*d_J); + timer.leave_subsection(); + + constraints.distribute(newton_update); + return std::make_pair(lin_it, lin_res); + } + template + void + Solid::update_end_timestep () + { + FilteredIterator::active_cell_iterator> + cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.begin_active()), + endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.end()); + for (; cell != endc; ++cell) + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + const std::vector > > lqph = + quadrature_point_history.get_data(cell); + Assert(lqph.size() == n_q_points, ExcInternalError()); + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) + lqph[q_point]->update_end_timestep(); + } + } + + template > + class FilteredDataOut : public DataOut + { + public: + FilteredDataOut (const unsigned int subdomain_id) + : + subdomain_id (subdomain_id) + {} + + virtual ~FilteredDataOut() {} + + virtual typename DataOut::cell_iterator + first_cell () + { + typename DataOut::active_cell_iterator + cell = this->dofs->begin_active(); + while ((cell != this->dofs->end()) && + (cell->subdomain_id() != subdomain_id)) + ++cell; + return cell; + } + + virtual typename DataOut::cell_iterator + next_cell (const typename DataOut::cell_iterator &old_cell) + { + if (old_cell != this->dofs->end()) + { + const IteratorFilters::SubdomainEqualTo predicate(subdomain_id); + return + ++(FilteredIterator + ::active_cell_iterator> + (predicate,old_cell)); + } + else + return old_cell; + } + + private: + const unsigned int subdomain_id; + }; + + template + void Solid::output_results(const unsigned int timestep, + const double current_time) const + { + // Output -> Ghosted vector + LA::MPI::BlockVector solution_total (locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator, + /*vector_writable = */ false); + LA::MPI::BlockVector residual (locally_owned_partitioning, + locally_relevant_partitioning, + mpi_communicator, + /*vector_writable = */ false); + solution_total = solution_n; + residual = system_rhs; + residual *= -1.0; + + // --- Additional data --- + Vector material_id; + Vector polynomial_order; + material_id.reinit(triangulation.n_active_cells()); + polynomial_order.reinit(triangulation.n_active_cells()); + std::vector partition_int (triangulation.n_active_cells()); + + FilteredDataOut data_out(this_mpi_process); + std::vector + data_component_interpretation(dim, + DataComponentInterpretation::component_is_part_of_vector); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + data_component_interpretation.push_back(DataComponentInterpretation::component_is_scalar); + + GridTools::get_subdomain_association (triangulation, partition_int); + + // Can't use filtered iterators here because the cell + // count "c" is incorrect for the parallel case + unsigned int c = 0; + typename DoFHandler::active_cell_iterator + cell = dof_handler.begin_active(), + endc = dof_handler.end(); + for (; cell!=endc; ++cell, ++c) + { + if (cell->subdomain_id() != this_mpi_process) continue; + + material_id(c) = static_cast(cell->material_id()); + } + + std::vector solution_name(n_components, "solution_"); + std::vector residual_name(n_components, "residual_"); + for (unsigned int c=0; c::type_dof_data, + data_component_interpretation); + data_out.add_data_vector(residual, + residual_name, + DataOut::type_dof_data, + data_component_interpretation); + const Vector partitioning(partition_int.begin(), + partition_int.end()); + data_out.add_data_vector (material_id, "material_id"); + data_out.add_data_vector (partitioning, "partitioning"); + data_out.build_patches(degree); + + struct Filename + { + static std::string get_filename_vtu (unsigned int process, + unsigned int timestep, + const unsigned int n_digits = 4) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << "." + << Utilities::int_to_string (process, n_digits) + << "." + << Utilities::int_to_string(timestep, n_digits) + << ".vtu"; + return filename_vtu.str(); + } + + static std::string get_filename_pvtu (unsigned int timestep, + const unsigned int n_digits = 4) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << "." + << Utilities::int_to_string(timestep, n_digits) + << ".pvtu"; + return filename_vtu.str(); + } + + static std::string get_filename_pvd (void) + { + std::ostringstream filename_vtu; + filename_vtu + << "solution-" + << (std::to_string(dim) + "d") + << ".pvd"; + return filename_vtu.str(); + } + }; + + // Write out main data file + const std::string filename_vtu = Filename::get_filename_vtu(this_mpi_process, timestep); + std::ofstream output(filename_vtu.c_str()); + data_out.write_vtu(output); + + // Collection of files written in parallel + // This next set of steps should only be performed + // by master process + if (this_mpi_process == 0) + { + // List of all files written out at this timestep by all processors + std::vector parallel_filenames_vtu; + for (unsigned int p=0; p < n_mpi_processes; ++p) + { + parallel_filenames_vtu.push_back(Filename::get_filename_vtu(p, timestep)); + } + + const std::string filename_pvtu (Filename::get_filename_pvtu(timestep)); + std::ofstream pvtu_master(filename_pvtu.c_str()); + data_out.write_pvtu_record(pvtu_master, + parallel_filenames_vtu); + + // Time dependent data master file + static std::vector > time_and_name_history; + time_and_name_history.push_back (std::make_pair (current_time, + filename_pvtu)); + const std::string filename_pvd (Filename::get_filename_pvd()); + std::ofstream pvd_output (filename_pvd.c_str()); + DataOutBase::write_pvd_record (pvd_output, time_and_name_history); + } + } + template + void Solid::compute_vertex_positions(std::vector &real_time, + std::vector > > &tracked_vertices, + const LA::MPI::BlockVector &solution_total) const + { + real_time.push_back(time.current()); + + std::vector vertex_found (tracked_vertices.size(), false); + std::vector > vertex_update (tracked_vertices.size()); + + FilteredIterator::active_cell_iterator> + cell (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.begin_active()), + endc (IteratorFilters::SubdomainEqualTo(this_mpi_process), + dof_handler.end()); + for (; cell != endc; ++cell) + { + Assert(cell->subdomain_id()==this_mpi_process, ExcInternalError()); + for (unsigned int v=0; v::vertices_per_cell; ++v) + { + for (unsigned int p=0; p pt_ref = tracked_vertices[p][0]; + if (cell->vertex(v).distance(pt_ref) < 1e-6*parameters.scale) + { + for (unsigned int d=0; dvertex_dof_index(v,u_block+d)); + + vertex_found[p] = true; + } + } + } + } + + for (unsigned int p=0; p0, ExcMessage("Vertex not found on any processor")); + Tensor<1,dim> update; + for (unsigned int d=0; d solid("parameters.prm"); + solid.run(); + } + catch (std::exception &exc) + { + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + 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 (...) + { + if (Utilities::MPI::this_mpi_process(MPI_COMM_WORLD) == 0) + { + std::cerr << std::endl << std::endl + << "----------------------------------------------------" + << std::endl; + std::cerr << "Unknown exception!" << std::endl << "Aborting!" + << std::endl + << "----------------------------------------------------" + << std::endl; + return 1; + } + } + return 0; +}