From 8f29ac02f917ceca4afbdd309b370e5c1fc50569 Mon Sep 17 00:00:00 2001 From: Wolfgang Bangerth Date: Fri, 14 Jan 2011 00:07:03 +0000 Subject: [PATCH] Add step-44 from Jean-Paul Pelteret and Andrew McBride. git-svn-id: https://svn.dealii.org/trunk@23187 0785d39b-7218-0410-832d-ea1e28bc413d --- deal.II/examples/step-44/Makefile | 150 ++ deal.II/examples/step-44/doc/builds-on | 1 + deal.II/examples/step-44/doc/intro.dox | 2 + deal.II/examples/step-44/doc/kind | 1 + deal.II/examples/step-44/doc/results.dox | 2 + deal.II/examples/step-44/doc/tooltip | 1 + deal.II/examples/step-44/step-44.cc | 2055 ++++++++++++++++++++++ 7 files changed, 2212 insertions(+) create mode 100644 deal.II/examples/step-44/Makefile create mode 100644 deal.II/examples/step-44/doc/builds-on create mode 100644 deal.II/examples/step-44/doc/intro.dox create mode 100644 deal.II/examples/step-44/doc/kind create mode 100644 deal.II/examples/step-44/doc/results.dox create mode 100644 deal.II/examples/step-44/doc/tooltip create mode 100644 deal.II/examples/step-44/step-44.cc diff --git a/deal.II/examples/step-44/Makefile b/deal.II/examples/step-44/Makefile new file mode 100644 index 0000000000..e54591f1af --- /dev/null +++ b/deal.II/examples/step-44/Makefile @@ -0,0 +1,150 @@ +# $Id$ + + +# For the small projects Makefile, you basically need to fill in only +# four fields. +# +# The first is the name of the application. It is assumed that the +# application name is the same as the base file name of the single C++ +# file from which the application is generated. +target = $(basename $(shell echo step-*.cc)) + +# The second field determines whether you want to run your program in +# debug or optimized mode. The latter is significantly faster, but no +# run-time checking of parameters and internal states is performed, so +# you should set this value to `on' while you develop your program, +# and to `off' when running production computations. +debug-mode = on + + +# As third field, we need to give the path to the top-level deal.II +# directory. You need to adjust this to your needs. Since this path is +# probably the most often needed one in the Makefile internals, it is +# designated by a single-character variable, since that can be +# reference using $D only, i.e. without the parentheses that are +# required for most other parameters, as e.g. in $(target). +D = ../../ + + +# The last field specifies the names of data and other files that +# shall be deleted when calling `make clean'. Object and backup files, +# executables and the like are removed anyway. Here, we give a list of +# files in the various output formats that deal.II supports. +clean-up-files = *gmv *gnuplot *gpl *eps *pov + + + + +# +# +# Usually, you will not need to change anything beyond this point. +# +# +# The next statement tell the `make' program where to find the +# deal.II top level directory and to include the file with the global +# settings +include $D/common/Make.global_options + + +# Since the whole project consists of only one file, we need not +# consider difficult dependencies. We only have to declare the +# libraries which we want to link to the object file, and there need +# to be two sets of libraries: one for the debug mode version of the +# application and one for the optimized mode. Here we have selected +# the versions for 2d. Note that the order in which the libraries are +# given here is important and that your applications won't link +# properly if they are given in another order. +# +# You may need to augment the lists of libraries when compiling your +# program for other dimensions, or when using third party libraries +libs.g = $(lib-deal2.g) +libs.o = $(lib-deal2.o) + + +# We now use the variable defined above which switch between debug and +# optimized mode to select the set of libraries to link with. Included +# in the list of libraries is the name of the object file which we +# will produce from the single C++ file. Note that by default we use +# the extension .g.o for object files compiled in debug mode and .o for +# object files in optimized mode (or whatever the local default on your +# system is instead of .o). +ifeq ($(debug-mode),on) + libraries = $(target).g.$(OBJEXT) $(libs.g) +else + libraries = $(target).$(OBJEXT) $(libs.o) +endif + + +# Now comes the first production rule: how to link the single object +# file produced from the single C++ file into the executable. Since +# this is the first rule in the Makefile, it is the one `make' selects +# if you call it without arguments. +$(target) : $(libraries) + @echo ============================ Linking $@ + @$(CXX) -o $@$(EXEEXT) $^ $(LIBS) $(LDFLAGS) + + +# To make running the application somewhat independent of the actual +# program name, we usually declare a rule `run' which simply runs the +# program. You can then run it by typing `make run'. This is also +# useful if you want to call the executable with arguments which do +# not change frequently. You may then want to add them to the +# following rule: +run: $(target) + @echo ============================ Running $< + @./$(target)$(EXEEXT) + + +# As a last rule to the `make' program, we define what to do when +# cleaning up a directory. This usually involves deleting object files +# and other automatically created files such as the executable itself, +# backup files, and data files. Since the latter are not usually quite +# diverse, you needed to declare them at the top of this file. +clean: + -rm -f *.$(OBJEXT) *~ Makefile.dep $(target)$(EXEEXT) $(clean-up-files) + + +# Since we have not yet stated how to make an object file from a C++ +# file, we should do so now. Since the many flags passed to the +# compiler are usually not of much interest, we suppress the actual +# command line using the `at' sign in the first column of the rules +# and write the string indicating what we do instead. +./%.g.$(OBJEXT) : + @echo ==============debug========= $( $@ \ + || (rm -f $@ ; false) + @if test -s $@ ; then : else rm $@ ; fi + + +# To make the dependencies known to `make', we finally have to include +# them: +include Makefile.dep + + diff --git a/deal.II/examples/step-44/doc/builds-on b/deal.II/examples/step-44/doc/builds-on new file mode 100644 index 0000000000..90d9d634fe --- /dev/null +++ b/deal.II/examples/step-44/doc/builds-on @@ -0,0 +1 @@ +step-18 diff --git a/deal.II/examples/step-44/doc/intro.dox b/deal.II/examples/step-44/doc/intro.dox new file mode 100644 index 0000000000..c59cfdf609 --- /dev/null +++ b/deal.II/examples/step-44/doc/intro.dox @@ -0,0 +1,2 @@ + +

Introduction

diff --git a/deal.II/examples/step-44/doc/kind b/deal.II/examples/step-44/doc/kind new file mode 100644 index 0000000000..56e049c91a --- /dev/null +++ b/deal.II/examples/step-44/doc/kind @@ -0,0 +1 @@ +solids diff --git a/deal.II/examples/step-44/doc/results.dox b/deal.II/examples/step-44/doc/results.dox new file mode 100644 index 0000000000..b5eaba9377 --- /dev/null +++ b/deal.II/examples/step-44/doc/results.dox @@ -0,0 +1,2 @@ +

Results

+ diff --git a/deal.II/examples/step-44/doc/tooltip b/deal.II/examples/step-44/doc/tooltip new file mode 100644 index 0000000000..2e08a00fc9 --- /dev/null +++ b/deal.II/examples/step-44/doc/tooltip @@ -0,0 +1 @@ +Quasi-static finite-strain elasticity diff --git a/deal.II/examples/step-44/step-44.cc b/deal.II/examples/step-44/step-44.cc new file mode 100644 index 0000000000..081f0404bb --- /dev/null +++ b/deal.II/examples/step-44/step-44.cc @@ -0,0 +1,2055 @@ + +/* Authors: Jean-Paul Pelteret, University of Cape Town, */ +/* Andrew McBride, University of Erlangen-Nuremberg, 2010 */ +/* */ +/* Copyright (C) 2010 by the deal.II authors */ +/* & Jean-Paul Pelteret and Andrew McBride */ +/* */ +/* This file is subject to QPL and may not be distributed */ +/* without copyright and license information. Please refer */ +/* to the file deal.II/doc/license.html for the text and */ +/* further information on this license. */ + +#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 + +using namespace dealii; + + // @sect3{Run-time parameters} +namespace Parameters +{ + // Finite Element system + struct FESystem + { + int poly_degree; + 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", + "1", + Patterns::Integer(), + "Displacement system polynomial order"); + + prm.declare_entry("Quadrature order", + "2", + Patterns::Integer(), + "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(); + } + + // Geometry + struct Geometry + { + int global_refinement; + double scale; + double p_p0; + + static void declare_parameters (ParameterHandler &prm); + void parse_parameters (ParameterHandler &prm); + }; + + void Geometry::declare_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + prm.declare_entry("Global refinement", + "2", + Patterns::Integer(), + "Global refinement level"); + + prm.declare_entry("Grid scale", + "1.0", + Patterns::Double(), + "Global grid scaling factor"); + + prm.declare_entry("Pressure ratio p/p0", + "40", + Patterns::Selection("20|40|60|80|100"), + "Ratio of applied pressure to reference pressure"); + } + prm.leave_subsection(); + } + + void Geometry::parse_parameters (ParameterHandler &prm) + { + prm.enter_subsection("Geometry"); + { + global_refinement = prm.get_integer("Global refinement"); + scale = prm.get_double("Grid scale"); + p_p0= prm.get_double("Pressure ratio p/p0"); + } + prm.leave_subsection(); + } + + // Materials + 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.49", + Patterns::Double(), + "Poisson's ratio"); + + prm.declare_entry("Shear modulus", + "1.0e6", + 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(); + } + + // Linear solver + struct LinearSolver + { + std::string type_lin; + double tol_lin; + double max_iterations_lin; + double ssor_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(), + "Linear solver residual (scaled by residual norm)"); + + prm.declare_entry("Max iteration multiplier", + "2", + Patterns::Double(), + "Linear solver iterations (multiples of the system matrix size)"); + + prm.declare_entry("SSOR Relaxation", + "0.6", + Patterns::Double(), + "SSOR 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"); + ssor_relaxation = prm.get_double("SSOR Relaxation"); + } + prm.leave_subsection(); + } + + // 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(), + "Number of Newton-Raphson iterations allowed"); + + prm.declare_entry("Tolerance force", + "1.0e-9", + Patterns::Double(), + "Force residual tolerance"); + + prm.declare_entry("Tolerance displacement", + "1.0e-3", + Patterns::Double(), + "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(); + } + + // Time + struct Time + { + double end_time; + double delta_t; + + 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(); + } + + // All parameters + 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{General tools} +namespace AdditionalTools +{ + template + SymmetricTensor<4,dim> outer_product_T23 (const SymmetricTensor<2,dim> & A, + const SymmetricTensor<2,dim> & B) + { + SymmetricTensor<4,dim> A_ik_B_jl; + + for (unsigned int i=0; i + void extract_submatrix(const std::vector< unsigned int > &row_index_set, + const std::vector< unsigned int > &column_index_set, + const MatrixType &matrix, + FullMatrix< double > &sub_matrix ) + { + + const unsigned int n_rows_submatrix = row_index_set.size(); + const unsigned int n_cols_submatrix = column_index_set.size(); + + sub_matrix.reinit(n_rows_submatrix, n_cols_submatrix); + + for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { + const unsigned int row = row_index_set[sub_row]; + Assert (row<=matrix.m(), ExcInternalError()); + + for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) { + const unsigned int col = column_index_set[sub_col]; + Assert (col<=matrix.n(), ExcInternalError()); + + sub_matrix(sub_row,sub_col) = matrix(row, col); + } + } + } + + template + void replace_submatrix(const std::vector< unsigned int > &row_index_set, + const std::vector< unsigned int > &column_index_set, + const MatrixType &sub_matrix, + FullMatrix< double > &matrix) + { + const unsigned int n_rows_submatrix = row_index_set.size(); + Assert (n_rows_submatrix<=sub_matrix.m(), ExcInternalError()); + const unsigned int n_cols_submatrix = column_index_set.size(); + Assert (n_cols_submatrix<=sub_matrix.n(), ExcInternalError()); + + for (unsigned int sub_row = 0; sub_row < n_rows_submatrix; ++sub_row) { + const unsigned int row = row_index_set[sub_row]; + Assert (row<=matrix.m(), ExcInternalError()); + + for (unsigned int sub_col = 0; sub_col < n_cols_submatrix; ++sub_col) { + const unsigned int col = column_index_set[sub_col]; + Assert (col<=matrix.n(), ExcInternalError()); + + matrix(row, col) = sub_matrix(sub_row, sub_col); + + } + } + } + +} + + // @sect3{Time class} +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 (void) {} + + const double & current (void) const {return time_current;} + const double & end (void) const {return time_end;} + const double & get_delta_t (void) const {return delta_t;} + const unsigned int & get_timestep (void) const {return timestep;} + void increment (void) {time_current += delta_t; ++timestep;} + +private: + unsigned int timestep; + double time_current; + const double time_end; + const double delta_t; +}; + + // @sect3{Neo-Hookean material} +template + class Material_NH +{ +public: + /// \brief Class constructor + Material_NH (const double & lambda, + const double & mu) + : + lambda_0 (lambda), + mu_0 (mu), + kappa_0 (lambda + 2.0/3.0*mu) + { } + virtual ~Material_NH (void) {}; + + // Stress and constitutive tensors + virtual SymmetricTensor<2, dim> get_T (const double & J, + const SymmetricTensor <2, dim> & B) + { + const double dW_dJ = get_dU_dtheta (J); + return mu_0*B + dW_dJ*J*I; + } + + virtual SymmetricTensor<4, dim> get_JC (const double & J, + const SymmetricTensor <2, dim> & B) + { + const double dW_dJ = get_dU_dtheta (J); + const double d2W_dJ2 = get_d2U_dtheta2 (J); + return J*( (dW_dJ + J*d2W_dJ2)*IxI - (2.0*dW_dJ)*II ); + } + + + // Volumetric quantities methods + double get_dU_dtheta (const double & d) {return kappa_0*(d - 1.0/d);} + double get_d2U_dtheta2 (const double & d) {return kappa_0*(1.0 + 1.0/(d*d));} + +protected: + // Material properties + const double mu_0; // Shear modulus + const double lambda_0; // Lame modulus + const double kappa_0; // Bulk modulus + + static SymmetricTensor<2, dim> const I; + static SymmetricTensor<4, dim> const IxI; + static SymmetricTensor<4, dim> const II; +}; + +template SymmetricTensor<2, dim> const Material_NH::I = SymmetricTensor<2, dim> (unit_symmetric_tensor ()); +template SymmetricTensor<4, dim> const Material_NH::IxI = SymmetricTensor<4, dim> (outer_product (I, I)); +template SymmetricTensor<4, dim> const Material_NH::II = SymmetricTensor<4, dim> (identity_tensor ()); + + // @sect3{Quadrature point history} +template + class PointHistory +{ +public: + PointHistory (void) + : + material (NULL), + pressure_n (0.0), + dilatation_n (1.0) + { } + virtual ~PointHistory (void) {delete material;} + + void setup_lqp ( Parameters::AllParameters & parameters ) + { + const double lambda = 2.0*parameters.mu*parameters.nu / (1.0-2.0*parameters.nu); + material = new Material_NH (lambda, + parameters.mu); + + // Initialise all tensors correctly + update_values (Tensor <2,dim> (), 0.0, 1.0); + } + + // Total Variables + void update_values (const Tensor<2, dim> & grad_u_n, + const double & pressure, + const double & dilatation) + { + // Calculated variables from displacement, displacement gradients + const Tensor <2,dim> F = static_cast > (unit_symmetric_tensor ()) + grad_u_n; + J = determinant(F); + F_inv = invert(F); + B_bar = std::pow(get_J(), -2.0/3.0) * symmetrize ( F* transpose (F) ); + + // Precalculated pressure, dilatation + pressure_n = pressure; + dilatation_n = dilatation; + + // Now that all the necessary variables are set, we can update the stress tensors + // Stress update can only update the stresses once the + // dilatation has been set as p = p(d) + T_bar = material->get_T (get_J(), get_B_bar()); + T_iso = dev_P*get_T_bar(); // Note: T_iso depends on T_bar + T_vol = get_pressure()*get_J()*I; + } + + // Displacement and strain + const double & get_dilatation(void) const {return dilatation_n;} + const double & get_J (void) const {return J;} + const Tensor <2,dim> & get_F_inv (void) const {return F_inv;} + const SymmetricTensor <2,dim> & get_B_bar (void) const {return B_bar;} + + // Volumetric terms + double get_dU_dtheta (void) { + return material->get_dU_dtheta(get_dilatation()); + } + + double get_d2U_dtheta2 (void) { + return material->get_d2U_dtheta2(get_dilatation()); + } + + // Stress + double get_pressure(void) {return pressure_n;} + const SymmetricTensor<2, dim> & get_T_iso (void) const {return T_iso;} + const SymmetricTensor<2, dim> & get_T_vol (void) const {return T_vol;}; + + // Tangent matrices + SymmetricTensor <4,dim> get_C_iso(void) + { + const double & J = get_J(); + const SymmetricTensor<2, dim> & B_bar = get_B_bar(); + const SymmetricTensor<2, dim> & T_iso = get_T_iso(); + + const SymmetricTensor <4,dim> T_iso_x_I = outer_product(T_iso, I); + const SymmetricTensor <4,dim> I_x_T_iso = outer_product(I, T_iso); + const SymmetricTensor <4,dim> CC_bar = material->get_JC (J, B_bar); + + return 2.0/3.0*trace(get_T_bar())*dev_P + - 2.0/3.0*(T_iso_x_I + I_x_T_iso) + + dev_P*CC_bar*dev_P; + } + + SymmetricTensor <4,dim> get_C_vol(void) + { + const double & p = get_pressure(); + const double & J = get_J(); + return p*J*(IxI - 2.0*II); + } + +private: + // === MATERIAL === + Material_NH * material; + + // ==== VOLUME, DISPLACEMENT AND STRAIN VARIABLES ==== + double dilatation_n; // Current dilatation + double J; + Tensor <2,dim> F_inv; + SymmetricTensor <2,dim> B_bar; + SymmetricTensor <2,dim> E; + + // ==== STRESS VARIABLES ==== + double pressure_n; // Current pressure + SymmetricTensor<2, dim> T_bar; + SymmetricTensor<2, dim> T_iso; + SymmetricTensor<2, dim> T_vol; + const SymmetricTensor<2, dim> & get_T_bar (void) const {return T_bar;} + + // Basis tensors + static SymmetricTensor<2, dim> const I; + static SymmetricTensor<4, dim> const IxI; + static SymmetricTensor<4, dim> const II; + static SymmetricTensor<4, dim> const dev_P; +}; + +template SymmetricTensor<2,dim> const PointHistory::I + = SymmetricTensor<2,dim> (unit_symmetric_tensor ()); +template SymmetricTensor<4,dim> const PointHistory::IxI + = SymmetricTensor<4,dim> (outer_product (I, I)); +template SymmetricTensor<4,dim> const PointHistory::II + = SymmetricTensor<4,dim> (identity_tensor ()); +template SymmetricTensor<4,dim> const PointHistory::dev_P + = SymmetricTensor<4,dim> (II - 1.0/3.0*IxI); + + + // @sect3{Quasi-static quasi-incompressible finite-strain solid} +template + class Solid +{ +public: + Solid (const std::string & input_file); + virtual ~Solid (void); + void run (void); + +private: + + // === DATA STRUCTS === + + struct 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 (void) { + cell_matrix = 0.0; + } + }; + + struct ScratchData_K + { + FEValues fe_values_ref; + + std::vector < std::vector< double > > Nx; + std::vector < std::vector< Tensor<2, dim> > > grad_Nx; + std::vector < std::vector< SymmetricTensor<2, dim> > > 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< double >(fe_cell.dofs_per_cell)), + grad_Nx (qf_cell.size(), + std::vector< Tensor<2, dim> >(fe_cell.dofs_per_cell)), + symm_grad_Nx (qf_cell.size(), + std::vector< SymmetricTensor<2, dim> >(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 (void) { + for (unsigned int q_point=0; q_point < grad_Nx.size(); ++q_point) { + for (unsigned int k=0; k < Nx.size(); ++k) { + Nx[q_point][k] = 0.0; + grad_Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + } + } + } + + }; + + struct PerTaskData_F + { + Vector cell_rhs; + std::vector local_dof_indices; + + PerTaskData_F (const unsigned int dofs_per_cell) + : + cell_rhs (dofs_per_cell), + local_dof_indices (dofs_per_cell) + { } + + void reset (void) { cell_rhs = 0.0; } + }; + + struct ScratchData_F + { + FEValues fe_values_ref; + FEFaceValues fe_face_values_ref; + + std::vector < std::vector< double > > Nx; + std::vector < std::vector< SymmetricTensor<2, dim> > > symm_grad_Nx; + std::vector< Vector > rhs_values; + + // Solution data + std::vector< std::vector > > solution_grads; + + ScratchData_F ( 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< double >(fe_cell.dofs_per_cell)), + symm_grad_Nx (qf_cell.size(), + std::vector< SymmetricTensor<2, dim> >(fe_cell.dofs_per_cell)), + rhs_values (qf_cell.size(), + Vector(dim)) + { } + + ScratchData_F ( const ScratchData_F & 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), + rhs_values (rhs.rhs_values) + { } + + void reset (void) { + for (unsigned int q_point=0; q_point < symm_grad_Nx.size(); ++q_point) { + for (unsigned int k=0; k < symm_grad_Nx[q_point].size(); ++k) { + Nx[q_point][k] = 0.0; + symm_grad_Nx[q_point][k] = 0.0; + rhs_values[q_point] = 0.0; + } + } + } + + }; + + struct PerTaskData_SC + { + FullMatrix cell_matrix; + std::vector local_dof_indices; + + // Calculation matrices (auto resized) + FullMatrix K_orig; + FullMatrix K_pu; + FullMatrix K_pt; + FullMatrix K_tt; + // Calculation matrices (manual resized) + FullMatrix K_pt_inv; + FullMatrix K_tt_inv; + FullMatrix K_con; + FullMatrix A; + FullMatrix B; + FullMatrix C; + + PerTaskData_SC (const unsigned int & dofs_per_cell, + const unsigned int & n_u, + const unsigned int & n_p, + const unsigned int & n_t) + : + cell_matrix (dofs_per_cell, + dofs_per_cell), + local_dof_indices (dofs_per_cell), + K_pt_inv (n_t, n_p), + K_tt_inv (n_t, n_t), + K_con (n_u, n_u), + A (n_t, n_u), + B (n_t, n_u), + C (n_p, n_u) + { } + + // Choose not to reset any data + // The matrix extraction and replacement tools will take care of this + void reset(void) { } + }; + + // Dummy struct for TBB + struct ScratchData_SC + { + ScratchData_SC (void) { } + ScratchData_SC (const ScratchData_SC & rhs) { } + void reset (void) { } + }; + + // Dummy struct for TBB + struct PerTaskData_UQPH + { + PerTaskData_UQPH (void) { } + void reset(void) { } + }; + + struct ScratchData_UQPH + { + const BlockVector & solution_total; + + std::vector< Tensor< 2, dim> > solution_grads_u_total; + std::vector solution_values_p_total; + std::vector solution_values_t_total; + + FEValues fe_values_ref; + + ScratchData_UQPH (const FiniteElement & fe_cell, + const QGauss & qf_cell, + const UpdateFlags uf_cell, + const BlockVector & solution_total) + : + fe_values_ref (fe_cell, + qf_cell, + uf_cell), + solution_grads_u_total (qf_cell.size()), + solution_values_p_total (qf_cell.size()), + solution_values_t_total (qf_cell.size()), + solution_total (solution_total) + { } + + ScratchData_UQPH (const ScratchData_UQPH & rhs) + : + fe_values_ref (rhs.fe_values_ref.get_fe(), + rhs.fe_values_ref.get_quadrature(), + rhs.fe_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_t_total (rhs.solution_values_t_total), + solution_total (rhs.solution_total) + { } + + void reset (void) + { + // Is this necessary? Won't the call to fe_values.get_gradient overwrite this data? + for (unsigned int q=0; q < qf_cell.size(); ++q) + { + solution_grads_u_total[q] = 0.0; + solution_values_p_total[q] = 0.0; + solution_values_t_total[q] = 0.0; + } + } + }; + + // === METHODS === + + /// \brief Print out a greeting for the user + void make_grid (void); + /// \brief Setup the Finite Element system to be solved + void system_setup (void); + void determine_component_extractors(void); + + /// \brief Assemble the system and right hand side matrices using multi-threading + void assemble_system_K (void); + void assemble_system_K_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_F (void); + void assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_F & scratch, + PerTaskData_F & data); + void copy_local_to_global_F (const PerTaskData_F & data); + void assemble_SC (void); + void assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_SC & scratch, + PerTaskData_SC & data); + void copy_local_to_global_SC (const PerTaskData_SC & data); + /// \brief Apply Dirichlet boundary values + void make_constraints (const int & it_nr, + ConstraintMatrix & constraints); + + // /// \brief Setup the quadrature point history for each cell + void setup_qph(void); + // /// \brief Update the quadrature points stress and strain values, and fibre directions + 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) {} + /// \brief Solve for the displacement using a Newton-Rhapson method + void solve_nonlinear_timestep (BlockVector & solution_delta); + void solve_linear_system (BlockVector & newton_update); + + /// \brief Error measurement + double get_error_dil (void); + + // Solution + BlockVector get_solution_total (const BlockVector & solution_delta); + + // Postprocessing + void output_results(void); + + // === ATTRIBUTES === + // Parameters + Parameters::AllParameters parameters; + + // Geometry + Triangulation triangulation; // Describes the triangulation + + // Time + Time time; + TimerOutput timer; + + // === Quadrature points === + std::vector< PointHistory > quadrature_point_history; // Quadrature point history + + // === Finite element system === + const unsigned int degree; + const FESystem fe; // Describes the global FE system + DoFHandler dof_handler_ref; // Describes the degrees of freedom + + unsigned int dofs_per_cell; // Number of degrees of freedom on each cell + const FEValuesExtractors::Vector u_fe; + const FEValuesExtractors::Scalar p_fe; + const FEValuesExtractors::Scalar t_fe; + + // Block description + 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 t_component = dim + 1; + + enum {u_dof=0 , p_dof, t_dof}; + std::vector dofs_per_block; + std::vector element_indices_u; + std::vector element_indices_p; + std::vector element_indices_t; + + // === Quadrature === + QGauss qf_cell; // Cell quadrature formula + QGauss qf_face; // Face quadrature formula + unsigned int n_q_points; // Number of quadrature points in a cell + unsigned int n_q_points_f; // Number of quadrature points in a face + + // === Stiffness matrix setup ==== + ConstraintMatrix constraints; // Matrix to keep track of all constraints + BlockSparsityPattern sparsity_pattern; // Sparsity pattern for the stiffness matrix + BlockSparseMatrix system_matrix; // Global stiffness matrix + BlockVector system_rhs; // Holds the right hand side vector + BlockVector solution_n; // Holds the solution vector: Total displacement over all time-steps + BlockVector soln_error; // Holds the error vector +}; + + // @sect3{Implementation of the Solid class} + + // @sect4{Public interface} +template + Solid::Solid (const std::string & input_file) + : + parameters (input_file), + triangulation (Triangulation::maximum_smoothing), + time (parameters.end_time, parameters.delta_t), + dof_handler_ref (triangulation), + 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 + qf_cell (parameters.quad_order), + qf_face (parameters.quad_order), + dofs_per_block (n_blocks), + u_fe (first_u_component), + p_fe (p_component), + t_fe (t_component), + timer (std::cout, + TimerOutput::summary, + TimerOutput::wall_times) +{ + n_q_points = qf_cell.size(); + n_q_points_f = qf_face.size(); + dofs_per_cell = fe.dofs_per_cell; + determine_component_extractors(); +} + +template + Solid::~Solid (void) +{ + dof_handler_ref.clear (); +} + +template + void Solid::run (void) +{ + // Pre-processing + make_grid (); + system_setup (); + output_results (); // Output initial grid position + time.increment(); + + BlockVector solution_delta (dofs_per_block); + solution_delta.collect_sizes (); + + while (time.current() < time.end()) { + solution_delta = 0.0; + + // Solve step and update total solution vector + solve_nonlinear_timestep (solution_delta); + solution_n += solution_delta; + + output_results (); + time.increment(); + } +} + + // @sect4{Solid::make_grid} +template + void Solid::make_grid (void) +{ + GridGenerator::hyper_rectangle ( triangulation, + Point (0.0, 0.0, 0.0), + Point (1.0, 1.0, 1.0), + true ); + GridTools::scale (parameters.scale, triangulation); + + // Need to refine at least once for the indentation problem + if (parameters.global_refinement == 0) triangulation.refine_global (1); + else triangulation.refine_global (parameters.global_refinement); + + // Apply different BC's to a patch on the top surface + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + for (; cell!=endc; ++cell) + { + if (cell->at_boundary() == true) { + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) { + // Find faces on the +y surface + if ( cell->face(face)->at_boundary() == true + && cell->face(face)->center()[2] == 1.0*parameters.scale) + { + if ( cell->face(face)->center()[0] < 0.5*parameters.scale + && cell->face(face)->center()[1] < 0.5*parameters.scale) + { + cell->face(face)->set_boundary_indicator (6); // Set a new boundary id on a patch + } + } + } + } + } +} + + // @sect4{Solid::system_setup} +template + void Solid::system_setup (void) +{ + timer.enter_subsection ("Setup system"); + + // Number of components per block + std::vector block_component (n_components, u_dof); // Displacement + block_component[p_component] = p_dof; // Pressure + block_component[t_component] = t_dof; // Dilatation + + // Setup DOF handler + dof_handler_ref.distribute_dofs (fe); + DoFRenumbering::Cuthill_McKee (dof_handler_ref); + DoFRenumbering::component_wise (dof_handler_ref, block_component); + // Count number of dofs per block + 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 system matrix + system_matrix.clear (); + { + const unsigned int n_dofs_u = dofs_per_block[u_dof]; + const unsigned int n_dofs_p = dofs_per_block[p_dof]; + const unsigned int n_dofs_t = dofs_per_block[t_dof]; + + BlockCompressedSimpleSparsityPattern csp (n_blocks, n_blocks); + + csp.block(u_dof,u_dof).reinit (n_dofs_u, n_dofs_u); + csp.block(u_dof,p_dof).reinit (n_dofs_u, n_dofs_p); + csp.block(u_dof,t_dof).reinit (n_dofs_u, n_dofs_t); + + csp.block(p_dof,u_dof).reinit (n_dofs_p, n_dofs_u); + csp.block(p_dof,p_dof).reinit (n_dofs_p, n_dofs_p); + csp.block(p_dof,t_dof).reinit (n_dofs_p, n_dofs_t); + + csp.block(t_dof,u_dof).reinit (n_dofs_t, n_dofs_u); + csp.block(t_dof,p_dof).reinit (n_dofs_t, n_dofs_p); + csp.block(t_dof,t_dof).reinit (n_dofs_t, n_dofs_t); + csp.collect_sizes(); + + DoFTools::make_sparsity_pattern (dof_handler_ref, csp, constraints, false); + sparsity_pattern.copy_from (csp); + } + + // the global system matrix will have the following structure + // | K'_uu | K_up | 0 | | dU_u | | dR_u | + // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | + // | 0 | K_tp | K_tt | | dU_t | | dR_t | + // reflect this structure in the sparsity pattern + 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) { + if( (ii < p_component) && (jj == t_component)) { + coupling[ii][jj] = DoFTools::none; + coupling[jj][ii] = DoFTools::none; + } else { + coupling[ii][jj] = DoFTools::always; + } + } + } + + system_matrix.reinit (sparsity_pattern); + + // Setup storage vectors + system_rhs.reinit (dofs_per_block); + system_rhs.collect_sizes (); + + solution_n.reinit (dofs_per_block); + solution_n.collect_sizes (); + solution_n.block(t_dof) = 1.0; // Dilatation is 1 in the initial configuration + + soln_error.reinit (dofs_per_block); + soln_error.collect_sizes (); + + // Set up the quadrature point history + setup_qph (); + + timer.leave_subsection(); +} + +// A way to extract subblocks from the matrix +template + void Solid::determine_component_extractors(void) +{ + element_indices_u.clear(); + element_indices_p.clear(); + element_indices_t.clear(); + + for (unsigned int k=0; k < fe.dofs_per_cell; ++k) { + // 0 = u, 1 = p, 2 = dilatation interpolation fields + const unsigned int k_group = fe.system_to_base_index(k).first.first; + if (k_group == u_dof) { + element_indices_u.push_back(k); + } + else if (k_group == p_dof) { + element_indices_p.push_back(k); + } + else if (k_group == t_dof) { + element_indices_t.push_back(k); + } + else { + Assert (k_group <= t_dof, ExcInternalError()); + } + } +} + + // @sect4{Solid::setup_qph} +template + void Solid::setup_qph (void) +{ + std::cout << " Setting up quadrature point data..." << std::endl; + + { + typename Triangulation::active_cell_iterator + cell = triangulation.begin_active(), + endc = triangulation.end(); + + unsigned int our_cells = 0; + for (; cell != endc; ++cell) { + cell->clear_user_pointer(); + ++our_cells; + } + + { + std::vector > tmp; + tmp.swap(quadrature_point_history); + } + + quadrature_point_history.resize(our_cells * n_q_points); + + unsigned int history_index = 0; + for (cell = triangulation.begin_active(); cell != endc; ++cell) { + cell->set_user_pointer(&quadrature_point_history[history_index]); + history_index += n_q_points; + } + + Assert(history_index == quadrature_point_history.size(), ExcInternalError()); + } + + // Setup initial data + typename DoFHandler::active_cell_iterator + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); + for (; cell != endc; ++cell) { + PointHistory* lqph = reinterpret_cast*> (cell->user_pointer()); + Assert(lqph >= &quadrature_point_history.front(), ExcInternalError()); + Assert(lqph < &quadrature_point_history.back(), ExcInternalError()); + + // Setup any initial information at displacement gauss points + for (unsigned int q_point = 0; q_point < n_q_points; ++q_point) { + lqph[q_point].setup_lqp( parameters ); + } + } +} + + // @sect4{Solid::update_qph_incremental} +template + void Solid::update_qph_incremental (const BlockVector & solution_delta) +{ + timer.enter_subsection("Update QPH data"); + std::cout << "Update QPH data..."<< std::endl; + + // Get total solution as it stands at this update increment + const BlockVector solution_total = get_solution_total(solution_delta); + const UpdateFlags uf_UQPH ( update_values | update_gradients ); + PerTaskData_UQPH per_task_data_UQPH; + ScratchData_UQPH scratch_data_UQPH (fe, + qf_cell, + uf_UQPH, + solution_total); + + 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(); +} + +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()); + Assert(scratch.solution_values_p_total.size() == n_q_points, ExcInternalError()); + Assert(scratch.solution_values_t_total.size() == n_q_points, ExcInternalError()); + + // Find the values and gradients at quadrature points inside the current cell + scratch.fe_values_ref.reinit(cell); + 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[t_fe].get_function_values (scratch.solution_total,scratch. solution_values_t_total); + + // === UPDATE DATA AT EACH GAUSS POINT === + // Update displacement and deformation gradient at all quadrature points + 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], + scratch.solution_values_p_total[q_point], + scratch.solution_values_t_total[q_point]); + } +} + + // @sect4{Solid::solve_nonlinear_timestep} +template + void Solid::solve_nonlinear_timestep (BlockVector & solution_delta) +{ +// timer.enter_subsection("Nonlinear solver"); + std::cout + << "Timestep " << time.get_timestep() + << std::endl; + + BlockVector newton_update (dofs_per_block); + newton_update.collect_sizes (); + + double res_u = 0.0, res_f = 0.0; + double res_u_0 = 1.0, res_f_0 = 1.0; + for (unsigned int it_nr=0; it_nr < parameters.max_iterations_NR; ++ it_nr) + { + std::cout + << std::endl + << "Newton iteration: " << it_nr + << std::endl; + + system_matrix = 0.0; + system_rhs = 0.0; + + // Check residual + assemble_system_K (); // Assemble stiffness matrix + assemble_system_F (); // Assemble RHS + make_constraints (it_nr, constraints); // Make boundary conditions + constraints.condense (system_matrix, system_rhs); // Apply BC's + + solve_linear_system (newton_update); + constraints.distribute(newton_update); // Populate the constrained DOF's with their values + + // Definition of residual for Newton's method: + // Newton's method: f(x + dx) = f(x) + f'(x).dx + ..... = 0 + // so: 0 = R + K.dU + // Sparsematrix residual: Mx = b + // r = b-Mx + // i.e.: r = -R - K.dU ~ 0 + system_matrix.residual (soln_error, newton_update, system_rhs); + res_u = newton_update.block(u_dof).l2_norm(); + res_f = soln_error.block(u_dof).l2_norm(); + + // Residual scaling factors + if (it_nr == 0) res_f_0 = res_f; + if (it_nr == 0) res_u_0 = res_u; + std::cout + << "Nonlinear system error: " + << std::endl << std::scientific + << " Solution update \t ||dU||: " << newton_update.l2_norm() + << "\t ||dU_u||: " << newton_update.block(u_dof).l2_norm() + << "\t ||dU_p||: " << newton_update.block(p_dof).l2_norm() + << "\t ||dU_t||: " << newton_update.block(t_dof).l2_norm() + << std::endl; + std::cout << std::scientific + << " Residual \t ||dF||: " << soln_error.l2_norm() + << "\t ||dR_u||: " << soln_error.block(u_dof).l2_norm() + << "\t ||dR_p||: " << soln_error.block(p_dof).l2_norm() + << "\t ||dR_t||: " << soln_error.block(t_dof).l2_norm() + << std::endl; + + // Check for solution convergence + if ( it_nr > 0 + && res_u/res_u_0 <= parameters.tol_u + && res_f/res_f_0 <= parameters.tol_f) + { + std::cout + << std::endl + << "Solution for timestep " << time.get_timestep() + << " converged on Newton iteration " << it_nr << "." + << std::endl + << "Relative displacement error: " << res_u/res_u_0 + << "\t Relative force error: " << res_f/res_f_0 + << "\t Dilatation error: " << get_error_dil() + << std::endl << std::endl; + +// timer.leave_subsection(); + return; + } + + // Current solution state unacceptable. Update and continue iterating. + solution_delta += newton_update; // Update current solution + update_qph_incremental (solution_delta); // Update quadrature point information + } + + throw(ExcMessage("No convergence in nonlinear solver!")); +} + +template + double Solid::get_error_dil (void) +{ + double v_e = 0.0; // Volume in current configuration + double V_e = 0.0; // Volume in reference configuration + + static FEValues fe_values_ref (fe, qf_cell, update_JxW_values); + + typename DoFHandler::active_cell_iterator + cell = dof_handler_ref.begin_active(), + endc = dof_handler_ref.end(); + for (; cell != endc; ++cell) { + fe_values_ref.reinit (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) { + v_e += lqph[q_point].get_dilatation() * fe_values_ref.JxW(q_point); + V_e += fe_values_ref.JxW(q_point); + } + } + + return std::abs((v_e - V_e)/V_e); // Difference between initial and current volume +} + +// Solution (valid at any Newton step) +template + BlockVector Solid::get_solution_total (const BlockVector & solution_delta) +{ + BlockVector solution_total (solution_n); + solution_total += solution_delta; + + return solution_total; +} + + // @sect4{Solid::solve_linear_system} +template + void Solid::solve_linear_system (BlockVector & newton_update) +{ + std::cout << "Solve linear system..." << std::endl; + + BlockVector A (dofs_per_block); + BlockVector B (dofs_per_block); + A.collect_sizes (); + B.collect_sizes (); + + // | K'_uu | K_up | 0 | | dU_u | | dR_u | + // K = | K_pu | K_tt^-1 | K_pt^-1 | , dU = | dU_p | , dR = | dR_p | + // | 0 | K_tp | K_tt | | dU_t | | dR_t | + + // Solve for du + { + // Do the static condensation to make K'_uu, and put K_pt^{-1} + // in the K_pt block and K_tt^{-1} in the K_pp block + assemble_SC(); + + // K'uu du = Ru − Kup Ktp^-1 (Rt − Ktt Kpt^{-1} Rp) + system_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), system_rhs.block(p_dof)); + system_matrix.block(t_dof, t_dof).vmult (B.block(t_dof), A.block(t_dof)); + A.block(t_dof).equ(1.0, system_rhs.block(t_dof), -1.0, B.block(t_dof)); + system_matrix.block(p_dof, t_dof).Tvmult(A.block(p_dof), A.block(t_dof)); + system_matrix.block(u_dof, p_dof).vmult(A.block(u_dof), A.block(p_dof)); + system_rhs.block(u_dof) -= A.block(u_dof); + + timer.enter_subsection("Linear solver"); + if (parameters.type_lin == "CG") + { + const int solver_its = system_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 < Vector > GVM; + SolverCG < Vector > solver_CG (solver_control, GVM); + + // SSOR -> much better than Jacobi for symmetric systems + PreconditionSSOR > preconditioner; + preconditioner.initialize (system_matrix.block(u_dof, u_dof), parameters.ssor_relaxation); + + solver_CG.solve (system_matrix.block(u_dof, u_dof), + newton_update.block(u_dof), + system_rhs.block(u_dof), + preconditioner); + + std::cout + << "\t Iterations: " << solver_control.last_step() + << "\n\t Residual: " << solver_control.last_value() + << std::endl; + } + else if (parameters.type_lin == "Direct") + { + SparseDirectUMFPACK A_direct; + A_direct.initialize(system_matrix.block(u_dof, u_dof)); + A_direct.vmult (newton_update.block(u_dof), + system_rhs.block(u_dof)); + } + else throw (ExcMessage("Linear solver type not implemented")); + timer.leave_subsection(); + } + + timer.enter_subsection("Linear solver postprocessing"); + // Postprocess for dp + { + // dp = Ktp^{-1} ( Rt − Ktt Kpt^{-1} (Rp − Kpu du) ) + system_matrix.block(p_dof, u_dof).vmult (A.block(p_dof), newton_update.block(u_dof)); + B.block(p_dof).equ(1.0, system_rhs.block(p_dof), -1.0, A.block(p_dof)); + system_matrix.block(p_dof, t_dof).vmult(A.block(t_dof), B.block(p_dof)); + system_matrix.block(t_dof, t_dof).vmult(B.block(t_dof), A.block(t_dof)); + A.block(t_dof).equ (1.0, system_rhs.block(t_dof), -1.0, B.block(t_dof)); + system_matrix.block(p_dof, t_dof).Tvmult (newton_update.block(p_dof), A.block(t_dof)); + } + + // Postprocess for dt + { + // dt = Ktt^{-1} (Rt − Ktp dp) + system_matrix.block(t_dof, p_dof).vmult (A.block(t_dof), newton_update.block(p_dof)); + system_rhs.block(t_dof) -= A.block(t_dof); + system_matrix.block(p_dof, p_dof).vmult (newton_update.block(t_dof), system_rhs.block(t_dof)); + } + timer.leave_subsection(); +} + + // @sect4{Solid::assemble_system_K} +template + void Solid::assemble_system_K (void) +{ + timer.enter_subsection("Assemble system matrix"); + std::cout << "Assemble system matrix..."<< std::endl; + + system_matrix = 0.0; // Clear the matrix + + const UpdateFlags uf_cell ( update_values | update_gradients | update_JxW_values ); + + PerTaskData_K per_task_data (dofs_per_cell); // Initialise members of per_task_data to the correct sizes. + ScratchData_K scratch_data (fe, qf_cell, uf_cell); + + WorkStream::run ( dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_system_K_one_cell, + &Solid::copy_local_to_global_K, + scratch_data, + per_task_data); + + timer.leave_subsection(); +} + +template + void Solid::copy_local_to_global_K (const PerTaskData_K & data) +{ + // Add the local contribution to the system matrix + for (unsigned int i=0; i + void Solid::assemble_system_K_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_K & scratch, + PerTaskData_K & data) +{ + data.reset(); // Reset data in the PerTaskData_K storage unit + scratch.reset(); // Reset data in the Scratch storage unit + scratch.fe_values_ref.reinit (cell); + cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have + PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); + + // Set up cell shape function gradients + static const SymmetricTensor<2, dim> I = unit_symmetric_tensor (); + 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 if (k_group == p_dof) { + scratch.Nx[q_point][k] = scratch.fe_values_ref[p_fe].value(k, q_point); + } + else if (k_group == t_dof) { + scratch.Nx[q_point][k] = scratch.fe_values_ref[t_fe].value(k, q_point); + } + else { + Assert (k_group <= t_dof, ExcInternalError()); + } + } + } + + // Build cell stiffness matrix + // Global and local system matrices are symmetric + // => Take advantage of this: Build only the lower half of the local matrix + // Only assemble 1/2 of the K_uu, K_pp = 0, K_tt blocks and the whole K_pt, K_ut, K_up blocks + for (unsigned int q_point=0; q_point < n_q_points; ++q_point) { + const Tensor <2,dim> T = static_cast < Tensor<2, dim> > (lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol()); + const SymmetricTensor <4,dim> C = lqph[q_point].get_C_iso() + lqph[q_point].get_C_vol(); + const double C_v = lqph[q_point].get_d2U_dtheta2(); + const double J = lqph[q_point].get_J(); + + const std::vector & N = scratch.Nx[q_point]; + const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; + const std::vector< Tensor <2,dim> > & B = 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; + + // Only assemble the lower diagonal part of the local matrix + 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 == j_group) && (i_group == u_dof ) ) { + data.cell_matrix(i,j) + += ( symm_B[i] * C * symm_B[j] // Material stiffness + + ( component_i == component_j ? + B[i][component_i] * T * B[j][component_j] : + 0.0 ) // Geometric stiffness. Only add this along local diagonals + ) * JxW; // K_uu + } + else if ( (i_group == p_dof) && (j_group == u_dof) ) { + data.cell_matrix(i,j) += N[i]*J*(symm_B[j]*I)*JxW; // K_pu + } + else if ( (i_group == t_dof) && (j_group == p_dof) ) { + data.cell_matrix(i,j) -= N[i]*N[j]*JxW; // K_tp + } + else if ( (i_group == j_group) && (i_group == t_dof) ) { + data.cell_matrix(i,j) += N[i]*C_v*N[j]*JxW; // K_tt + } + else Assert ((i_group <= t_dof) && (j_group <= t_dof), ExcInternalError()); + } // END j LOOP + } // END i LOOP + + } // END q_point LOOP + + // Global and local system matrices are symmetric + // => Copy the upper half of the local matrix in the bottom half of the local matrix + for (unsigned int i=0; i + void Solid::assemble_system_F (void) +{ + timer.enter_subsection("Assemble system RHS"); + std::cout << "Assemble system RHS..."<< std::endl; + + system_rhs = 0.0; // Clear the vector + + const UpdateFlags uf_cell ( update_values | update_gradients | update_JxW_values ); + const UpdateFlags uf_face ( update_values | update_normal_vectors | update_JxW_values); + + PerTaskData_F per_task_data (dofs_per_cell); // Initialise members of per_task_data to the correct sizes. + ScratchData_F 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_F_one_cell, + &Solid::copy_local_to_global_F, + scratch_data, + per_task_data ); + + timer.leave_subsection(); +} + +template + void Solid::copy_local_to_global_F (const PerTaskData_F & data) +{ + // Add the local contribution to the system RHS vector + for (unsigned int i=0; i + void Solid::assemble_system_F_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_F & scratch, + PerTaskData_F & data) +{ + data.reset(); // Reset data in the PerTaskData_K storage unit + scratch.reset(); // Reset data in the ScratchData_F storage unit + scratch.fe_values_ref.reinit (cell); + cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have + PointHistory *lqph = reinterpret_cast*>(cell->user_pointer()); + + // Precompute some data + 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 T = lqph[q_point].get_T_iso() + lqph[q_point].get_T_vol(); + const double J = lqph[q_point].get_J(); + const double D = lqph[q_point].get_dilatation(); + const double p = lqph[q_point].get_pressure(); + const double p_star = lqph[q_point].get_dU_dtheta(); + + const std::vector< double > & N = scratch.Nx[q_point]; + const std::vector< SymmetricTensor <2,dim> > & symm_B = scratch.symm_grad_Nx[q_point]; + const double JxW = scratch.fe_values_ref.JxW(q_point); + + for (unsigned int i=0; iat_boundary() == true) + { + static const Tensor <2, dim> I = static_cast < Tensor <2, dim> > ( unit_symmetric_tensor () ); + + for (unsigned int face=0; face < GeometryInfo::faces_per_cell; ++face) + { + if ( cell->face(face)->at_boundary() == true + && cell->face(face)->boundary_indicator() == 6 ) + { + 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); + + // Traction in reference configuration + // t_0 = p*N + static const double p0 = -4.0/(parameters.scale*parameters.scale); // Reference pressure of 4 Pa + const double time_ramp = (time.current() / time.end()); // Linearly ramp up the pressure with time + const double pressure = p0 * parameters.p_p0 * time_ramp; + const Tensor <1,dim> traction = pressure * N; + + for (unsigned int i=0; i < dofs_per_cell; ++i) { + // Determine the dimensional component that matches the dof component (i.e. i % dim) + 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); + + // Add traction vector contribution to the local RHS vector (displacement dofs only) + data.cell_rhs(i) += (Ni * traction[component_i]) // Contribution from external forces + * JxW; + } + } // END i LOOP + } // END face q_point LOOP + } // END at boundary check LOOP + + } // END face LOOP + } +} + + // @sect4{Solid::assemble_system_SC} +template + void Solid::assemble_SC (void) +{ + timer.enter_subsection("Perform static condensation"); + + PerTaskData_SC per_task_data (dofs_per_cell, + element_indices_u.size(), + element_indices_p.size(), + element_indices_t.size()); // Initialise members of per_task_data to the correct sizes. + ScratchData_SC scratch_data; + + WorkStream::run ( dof_handler_ref.begin_active(), + dof_handler_ref.end(), + *this, + &Solid::assemble_SC_one_cell, + &Solid::copy_local_to_global_SC, + scratch_data, + per_task_data ); + + timer.leave_subsection(); +} + +template + void Solid::copy_local_to_global_SC (const PerTaskData_SC & data) +{ + // Add the local contribution to the system matrix + for (unsigned int i=0; i + void Solid::assemble_SC_one_cell (const typename DoFHandler::active_cell_iterator & cell, + ScratchData_SC & scratch, + PerTaskData_SC & data) +{ + data.reset(); + scratch.reset(); + cell->get_dof_indices (data.local_dof_indices); // Find out which global numbers the degrees of freedom on this cell have + + // The local stifness matrix K_e is: + // | K_uu | K_up | 0 | + // | K_pu | 0 | K_pt | + // | 0 | K_tp | K_tt | + // + // We are going to exploit the zeros for post-processing as: + // | K'_uu | K_up | 0 | + // | K_pu | K_tt^-1 | K_pt^-1 | + // | 0 | K_tp | K_tt | + // with K'_uu = K_uu + Kup Ktp^{-1} Ktt Kpt^{-1} Kpu + + // NOTE: + // GLOBAL Data already exists in the K_uu, K_pt, K_tp subblocks + // + // For the K_uu block in particular, this means that contributions have been + // added from the surrounding cells, so we need to be careful when we manipulate this block. + // We can't just erase the subblocks and + // Additionally the copy_local_to_global operation is a "+=" operation -> need to take this + // into account + // + // So the intermediate matrix that we need to get from what we have in K_uu and what we + // are actually wanting is: + // | K'_uu - K_uu | 0 | 0 | + // | 0 | K_tt^-1 | K_pt^-1 - K_pt | + // | 0 | 0 | 0 | + // + // Strategy to get the subblocks we want: + // K'_uu: Since we don't have access to K_uu^h, but we know its contribution is added to the global + // K_uu matrix, we just want to add the element wise static-condensation + // K'_uu^h = K_uu^h + K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h + // Since we already have K_uu^h in the system matrix, we just need to do the following + // K'_uu^h == (K_uu^h += K_up^h K_tp^{-1}^h K_tt^h K_pt^{-1}^h K_pu^h) + // K_pt^-1: Similarly, K_pt exists in the subblock. Since the copy operation is a += operation, we need + // to subtract the existing K_pt submatrix in addition to "adding" that which we wish to + // replace it with. + // K_tp^-1: Same as above + // K_tt^-1: Nothing exists in the original K_pp subblock, so we can just add this contribution as is. + + // Extract element data from the system matrix + AdditionalTools::extract_submatrix(data.local_dof_indices, + data.local_dof_indices, + system_matrix, + data.K_orig); + AdditionalTools::extract_submatrix(element_indices_p, + element_indices_u, + data.K_orig, + data.K_pu); + AdditionalTools::extract_submatrix(element_indices_p, + element_indices_t, + data.K_orig, + data.K_pt); + AdditionalTools::extract_submatrix(element_indices_t, + element_indices_t, + data.K_orig, + data.K_tt); + + // Place K_pt^-1 in the K_pt block + data.K_pt_inv.invert(data.K_pt); + data.K_pt_inv.add (-1.0, data.K_pt); + AdditionalTools::replace_submatrix(element_indices_p, + element_indices_t, + data.K_pt_inv, + data.cell_matrix); + + // Place K_tt^-1 in the K_pp block + data.K_tt_inv.invert(data.K_tt); + AdditionalTools::replace_submatrix(element_indices_p, + element_indices_p, + data.K_tt_inv, + data.cell_matrix); + + // Make condensation terms to add to the K_uu block + data.K_pt_inv.mmult(data.A, data.K_pu); + data.K_tt.mmult(data.B, data.A); + data.K_pt_inv.Tmmult(data.C, data.B); // Symmetric matrix + data.K_pu.Tmmult(data.K_con, data.C); // Symmetric matrix + AdditionalTools::replace_submatrix(element_indices_u, + element_indices_u, + data.K_con, + data.cell_matrix); +} + + // @sect4{Solid::make_constraints} +template + void Solid::make_constraints (const int & it_nr, + ConstraintMatrix & constraints) +{ + std::cout << "Make constraints..."<< std::endl; + + constraints.clear(); + const bool apply_dirichlet_bc = (it_nr == 0); + + // Boundary conditions: + // b_id 0: -x face: Zero x-component of displacement : Symmetry plane + // b_id 2: -y face: Zero y-component of displacement : Symmetry plane + // b_id 4: -z face: Zero z-component of displacement : Symmetry plane + + // b_id 5: +z face: Zero x-component and Zero y-component + // b_id 6: Applied pressure face: Zero x-component and Zero y-component + // b_id 1: +x face: Traction free + // b_id 3: +y face: Traction free + { + const int boundary_id = 0; + + std::vector< bool > components (n_components, false); + components[0] = true; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + } + { + const int boundary_id = 2; + + std::vector< bool > components (n_components, false); + components[1] = true; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + } + { + const int boundary_id = 4; + std::vector< bool > components (n_components, false); + components[2] = true; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + } + { + const int boundary_id = 5; + std::vector< bool > components (n_components, true); + components[2] = false; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + } + { + const int boundary_id = 6; + std::vector< bool > components (n_components, true); + components[2] = false; + + if (apply_dirichlet_bc == true) { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + else { + VectorTools::interpolate_boundary_values ( dof_handler_ref, boundary_id, ZeroFunction(n_components), constraints, components ); + } + } + + constraints.close(); +} + + // @sect4{Solid::output_results} +template + void Solid::output_results(void) +{ + DataOut data_out; + + 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); + + std::vector solution_name (dim, "displacement"); + solution_name.push_back ("pressure"); + solution_name.push_back ("dilatation"); + + data_out.attach_dof_handler (dof_handler_ref); + data_out.add_data_vector (solution_n, + solution_name, + DataOut::type_dof_data, data_component_interpretation); +// MappingQEulerian q_mapping (degree, solution_n.block(u_dof), dof_handler_ref); +// MappingQEulerian q_mapping (degree, solution_n, dof_handler_ref); + Vector soln; + soln.reinit(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} +int main () +{ + 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; +} -- 2.39.5