deal.II version GIT relicensing-6842-g793a97d2aa 2026-10-02 14:00:01+00:00
\(\newcommand{\dealvcentcolon}{\mathrel{\mathop{:}}}\) \(\newcommand{\dealcoloneq}{\dealvcentcolon\mathrel{\mkern-1.2mu}=}\) \(\newcommand{\jump}[1]{\left[\!\left[ #1 \right]\!\right]}\) \(\newcommand{\average}[1]{\left\{\!\left\{ #1 \right\}\!\right\}}\)
Loading...
Searching...
No Matches
The 'L-BFGS phasefield solver ' code gallery program

This program was contributed by Tao Jin <[email protected]>.
It comes without any warranty or support by its authors or the authors of deal.II.

This program is part of the deal.II code gallery and consists of the following files (click to inspect):

Pictures from this code gallery program

Annotated version of README.md

Phasefield L-BFGS monolithic solver

A monolithic solver based on the limited-memory BFGS (L-BFGS) method for phase-field fracture simulations.

Purpose

This repository provides the source code and the input files for the numerical examples used in the paper titled “A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement”. The L-BFGS monolithic solver has the following features:

  1. It uses the limited-memory BFGS (L-BFGS) method to overcome the non-convexity of the total energy functional of the phase-field fracture formulation.
  2. It uses the history variable (the maximum of the positive strain energy in history) approach to enforce the phase-field irreversibility.
  3. It adopts an adaptive mesh refinement technique to reduce computational cost.
  4. It works for both 2D and 3D phase-field fracture simulations.

Content

The repository contains the following content:

  1. the source code of the L-BFGS method for the phase-field monolithic solver.
  2. the input files for several 2D and 3D phase-field fracture simulations included in the aforementioned manuscript.

Latest update:

  1. (Sept. 10th, 2025) Add a compiler macro so that the code works for the older version of deal.ii due to the interface change of the function interpolate() in the SolutionTransfer class.
  2. (Sept. 1st, 2025) Add a new gradient-based line search method. Comparing with the previously implemented Strong Wolfe line search, the new line search method could reduce the wall-clock time by 30 to 50 percent.

Representative results

A series of widely adopted phase-field crack benchmark problems are included in this code. Here are some examples:

  1. Simple tension test (pre-refined mesh): Tension test, prerefine mesh Tension test, mesh 15
  2. Simple shear test (adpatively refined mesh): Shear test, mesh 14 Shear test, solution 14
  3. 3D torsion test (adpatively refined mesh): Torsion test, mesh 3 Torsion test, solution 3

Phase-field model adopted in this work

The phase-field fracture approach aims to minimize the following energy functional: eq. 1 where the approximated crack surface is defined as eq. 2 the strain energy density function is based on the isotropic linear elasticity and the additive decomposition eq. 3 and the phase-field degradation function adopts the following qudratic form eq. 4

Using the divergence theorem and the technique of integration by parts, the corresponding Euler-Lagrange equations are written as eq. 5 In the actual implementation, however, the governing equations of the cracked solid system are modified as eq. 6 where a viscosity regularization term is introduced to stabilize the numerical treatment. Even though this regularization term might not be necessary, we still include it in the formulation for flexibility. This term can always be turned off by setting the viscosity coefficient as zero. Among various approaches to enforce the phase-field irreversibility, the approach based on the history variable of the positive strain energy is adopted due to its relative simplicity, eq. 8

Main idea of the limited-memory BFGS (L-BFGS) solver

The classical BFGS method (a type of quasi-Newton method) involves the following approximate Hessian matrix update during each iteration: eq. 7

The problem of this update in the context of finite element simulations is that the second term and the third term both generate a fully dense matrix of n by n needs to be stored, where n represents the number of degrees of freedom.** The above Hessian matrix update is too restrictive even for a mid-size finite element problem due to the required memory for the storage of the fully dense matrix. This limitation motivated this work to introduce the limited-memory feature for the phase-field crack simulations. The limited-memory BFGS method implemented in this work follows the algorithm represented in Chapter 7.2 (page 176) of the following textbook by Nocedal J, Wright SJ. Numerical optimization (2nd edition). Springer New York, NY, 2006.

How to compile

The L-BFGS finite element procedure is implemented in deal.II (originally with version 9.4.0 and also works for 9.5.1, it is also tested with the develop branch as Sept. 10th, 2025). In order to use the code (main.cc) provided here, deal.II should be configured with MPI and at least with the interfaces to BLAS, LAPACK, Threading Building Blocks (TBB), and UMFPACK. For optional interfaces to other software packages, see https://www.dealii.org/developer/readme.html.

Once the deal.II library is compiled, for instance, to "~/dealii-dev/bin/", follow the steps listed below:

  1. cmake -DDEAL_II_DIR=~/dealii-dev/bin/ .
  2. make debug or make release
  3. make

How to run

  1. Go into one of the examples folders.
  2. For instance, to run a 2D test case: go into examples/simple_tension_test/
  3. Run via ./../../main 2
  4. For instance, to run a 3D test case: go into examples/3D_torsion/
  5. Run via ./../../main 3

How to expand this code

If you want to use the current code to solve new 2D or 3D phase-field crack problems, you need to do the following:

  1. Add a new mesh under the function void make_grid().
  2. Add the boundary conditions for your new mesh in the function void make_constraints(const unsigned int it_nr).
  3. Modify the text file timeDataFile for the load/time step sizes and the text file materialDataFile for the material properties.
  4. Modify the input file parameters.prm accordingly.

If you want to use a new phase-field degradation function (the current code uses the standard quadratic degradation function), you can modify the following functions double degradation_function(const double d), double degradation_function_derivative(const double d), and double degradation_function_2nd_order_derivative(const double d).

If you want to modify the phase-field model completely but still use the L-BFGS monolithic solver, you need to modify the calculations of the initial BFGS matrix

void assemble_system_B0_one_cell(
const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM & scratch,
PerTaskData_ASM & data) const;

and the residuals

void assemble_system_rhs_BFGS_one_cell(
const typename DoFHandler<dim>::active_cell_iterator &cell,
ScratchData_ASM_RHS_BFGS & scratch,
PerTaskData_ASM_RHS_BFGS & data) const;

How to cite this work:

Jin T, Li Z, Chen K. A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement. Int J Numer Methods Eng. 2024;e7572. doi: 10.1002/nme.7572

@Article{2024:jin.li.ea:novel,
author = {Jin, Tao and Li, Zhao and Chen, Kuiying},
title = {A novel phase-field monolithic scheme for brittle crack propagation based on the limited-memory BFGS method with adaptive mesh refinement},
journal = {International Journal for Numerical Methods in Engineering},
year = 2024,
volume = 125,
number = 22,
pages = {e7572},
month = nov,
issn = {1097-0207},
url = {https://onlinelibrary.wiley.com/doi/10.1002/nme.7572},
doi = {10.1002/nme.7572},
publisher = {Wiley}
}

Annotated version of SpectrumDecomposition.cc

  #include "SpectrumDecomposition.h" // use double quote
  #include <deal.II/base/tensor.h>
  #include <deal.II/base/symmetric_tensor.h>
  #include <deal.II/base/exceptions.h>
 
 
  #include <iostream>
 
  namespace usr_spectrum_decomposition
  {
  void print()
  {
  Assert(3 >= 2, ::StandardExceptions::ExcNotImplemented());
  ::Tensor<1, 3> Br_tilde;
  std::cout << Br_tilde << std::endl;
  std::cout << "Hello world!" << std::endl;
  }
 
  double positive_ramp_function(const double x)
  {
  return std::fmax(x, 0.0);
  }
 
  double negative_ramp_function(const double x)
  {
  return std::fmin(x, 0.0);
  }
 
  double heaviside_function(const double x)
  {
  if (std::fabs(x) < 1.0e-16)
  return 0.5;
 
  if (x > 0)
  return 1.0;
  else
  return 0.0;
  }
  }
 
*  *  *  struct InterferenceTaperTransform *  
Definition ALE_Finite_Strain_Plasticity.h:7883
StandardExceptions::ExcNotImplemented
static ::ExceptionBase & ExcNotImplemented()
Assert
#define Assert(cond, exc)
Definition exception_macros.h:543

Annotated version of SpectrumDecomposition.h

  #ifndef usrcodes_spectrum_decomposition_h
  #define usrcodes_spectrum_decomposition_h
  #include <deal.II/base/tensor.h>
  #include <deal.II/base/symmetric_tensor.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/lapack_full_matrix.h>
  #include <deal.II/base/patterns.h>
  #include <fstream>
  #include <iostream>
 
 
  namespace usr_spectrum_decomposition
  {
  using namespace dealii;
  void print();
 
  double positive_ramp_function(const double x);
 
  double negative_ramp_function(const double x);
 
  double heaviside_function(const double x);
 

templated function has to be defined in the header file perform a spectrum decomposition of a symmetric tensor input: a symmetric tensor (SymmetricTensor<2, matrix_dimension>) output: eigenvalues (Vector<double>) eigenvectors (std::vector<Tensor<1, dim>>)

  template <int dim>
  void spectrum_decomposition(SymmetricTensor<2, dim> const & symmetric_tensor,
  Vector<double> & myEigenvalues,
  std::vector<Tensor<1, dim>> & myEigenvectors)
  {
 
  const std::array< std::pair< double, Tensor< 1, dim > >, dim >
  myEigenSystem = eigenvectors(symmetric_tensor);
 
  for (int i = 0; i < dim; i++)
  {
  myEigenvalues[i] = myEigenSystem[i].first;
  myEigenvectors[i] = myEigenSystem[i].second;
  }
  }
 
  template <int dim>
  SymmetricTensor<2, dim> positive_tensor(Vector<double> const & eigenvalues,
  std::vector<Tensor<1, dim>> const & eigenvectors)
  {
  SymmetricTensor<2, dim> positive_part_tensor;
  positive_part_tensor = 0;
  for (int i = 0; i < dim; i++)
  positive_part_tensor += positive_ramp_function(eigenvalues[i])
  * symmetrize(outer_product(eigenvectors[i],
  eigenvectors[i]));
  return positive_part_tensor;
  }
 
  template <int dim>
  SymmetricTensor<2, dim> negative_tensor(Vector<double> const & eigenvalues,
  std::vector<Tensor<1, dim>> const & eigenvectors)
  {
  SymmetricTensor<2, dim> negative_part_tensor;
  negative_part_tensor = 0;
  for (int i = 0; i < dim; i++)
  negative_part_tensor += negative_ramp_function(eigenvalues[i])
  * symmetrize(outer_product(eigenvectors[i],
  eigenvectors[i]));
  return negative_part_tensor;
  }
 
  template <int dim>
  void positive_negative_projectors(Vector<double> const & eigenvalues,
  std::vector<Tensor<1, dim>> const & eigenvectors,
  SymmetricTensor<4, dim> & positive_projector,
  SymmetricTensor<4, dim> & negative_projector)
  {
  Assert(dim <= 3,
  ExcMessage("Project tensors only work for dim <= 3."));
 
  std::array<SymmetricTensor<2, dim>, dim> M;
  for (int a = 0; a < dim; a++)
  M[a] = symmetrize(outer_product(eigenvectors[a], eigenvectors[a]));
 
  std::array<SymmetricTensor<4, dim>, dim> Q;
  for (int a = 0; a < dim; a++)
  Q[a] = outer_product(M[a], M[a]);
 
  std::array<std::array<SymmetricTensor<4, dim>, dim>, dim> G;
  for (int a = 0; a < dim; a++)
  for (int b = 0; b < dim; b++)
  for (int i = 0; i < dim; i++)
  for (int j = 0; j < dim; j++)
  for (int k = 0; k < dim; k++)
  for (int l = 0; l < dim; l++)
  G[a][b][i][j][k][l] = M[a][i][k] * M[b][j][l]
  + M[a][i][l] * M[b][j][k];
 
  positive_projector = 0;
  for (int a = 0; a < dim; a++)
  {
  double lambda_a = eigenvalues[a];
  positive_projector += heaviside_function(lambda_a)
  * Q[a];
  for (int b = 0; b < dim; b++)
  {
  if (b != a)
  {
  double lambda_b = eigenvalues[b];
  double v_ab = 0.0;
  if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
  v_ab = (positive_ramp_function(lambda_a) - positive_ramp_function(lambda_b))
  / (lambda_a - lambda_b);
  else
  v_ab = 0.5 * ( heaviside_function(lambda_a)
  + heaviside_function(lambda_b) );
  positive_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
  }
  }
  }
 
  negative_projector = 0;
  for (int a = 0; a < dim; a++)
  {
  double lambda_a = eigenvalues[a];
  negative_projector += heaviside_function(-lambda_a)
  * Q[a];
  for (int b = 0; b < dim; b++)
  {
  if (b != a)
  {
  double lambda_b = eigenvalues[b];
  double v_ab = 0.0;
  if (std::fabs(lambda_a - lambda_b) > 1.0e-12)
  v_ab = (negative_ramp_function(lambda_a) - negative_ramp_function(lambda_b))
  / (lambda_a - lambda_b);
  else
  v_ab = 0.5 * ( heaviside_function(-lambda_a)
  + heaviside_function(-lambda_b) );
  negative_projector += 0.5 * v_ab * 0.5 * (G[a][b] + G[b][a]);
  }
  }
  }
 
  }
 
  }
  #endif
Physics::Elasticity::Kinematics::l
Tensor< 2, dim, Number > l(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)
Physics::Elasticity::Kinematics::b
SymmetricTensor< 2, dim, Number > b(const Tensor< 2, dim, Number > &F)
symmetrize
constexpr SymmetricTensor< 2, dim, Number > symmetrize(const Tensor< 2, dim, Number > &t)
Definition symmetric_tensor.h:3572
eigenvalues
std::array< Number, 1 > eigenvalues(const SymmetricTensor< 2, 1, Number > &T)
positive_negative_projectors
std::tuple< SymmetricTensor< 2, dim, Number >, SymmetricTensor< 2, dim, Number >, SymmetricTensor< 4, dim, Number >, SymmetricTensor< 4, dim, Number > > positive_negative_projectors(const SymmetricTensor< 2, dim, Number > &original_tensor)
Definition symmetric_tensor.h:3448
eigenvectors
std::array< std::pair< Number, Tensor< 1, dim, Number > >, dim > eigenvectors(const SymmetricTensor< 2, dim, Number > &T, const SymmetricTensorEigenvectorMethod method=SymmetricTensorEigenvectorMethod::ql_implicit_shifts)
outer_product
constexpr SymmetricTensor< 4, dim, Number > outer_product(const SymmetricTensor< 2, dim, Number > &t1, const SymmetricTensor< 2, dim, Number > &t2)
Definition symmetric_tensor.h:3339

Annotated version of Utilities.h

  #ifndef usrcodes_utilities_h
  #define usrcodes_utilities_h
  #include <deal.II/grid/tria.h>
  #include <deal.II/dofs/dof_handler.h>
 
  #include <fstream>
  #include <iostream>
 
 
  namespace usr_utilities
  {
  using namespace dealii;
 
  template <int dim>
  std::vector<types::global_dof_index> get_vertex_dofs(
  const typename Triangulation<dim>::active_vertex_iterator &vertex,
  const DoFHandler<dim> &dof_handler)
  {
  DoFAccessor<0, dim, dim, false> vertex_dofs(
  &(dof_handler.get_triangulation()),
  vertex->level(),
  vertex->index(),
  &dof_handler);
  const unsigned int n_dofs = dof_handler.get_fe().dofs_per_vertex;
  std::vector<types::global_dof_index> dofs(n_dofs);
  for (unsigned int i = 0; i < n_dofs; ++i)
  {
  dofs[i] = vertex_dofs.vertex_dof_index(0, i);
  }
  return dofs;
  }
  }
  #endif

Annotated version of examples/L_shape_bending_cyclic_loading/plot_reaction_force.py

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# import matplotlib.ticker as mtick
time1, forceX1, forceY1 = np.loadtxt('Reaction_force.hist',
delimiter='\t', unpack=True)
#time2, forceX2, forceY2 = np.loadtxt('Reaction_force_timestep_10-5.hist',
# delimiter='\t', unpack=True)
font = {'size': 18}
matplotlib.rc('font', **font)
labels = []
labels.append("Loading")
labels.append("Unloading")
symbols = ['-k', '--*m', '-.b^', ':gx', '-rD', 'bx']
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:51], forceY1[:51], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(time1[51:], forceY1[51:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8,markevery=2)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Pseudo time')
plt.ylabel('Reaction force (kN)', multialignment='center')
#plt.xlim(0,0.05)
# plt.xticks(np.arange(0, 6, step=1))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend()
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_time.eps", bbox_inches='tight', pad_inches=0.1)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:51], forceY1[:51], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(2.0 - time1[51:], forceY1[51:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8,markevery=2)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Displacement (mm)')
plt.ylabel('Reaction force (kN)', multialignment='center')
#plt.xlim(0,0.05)
# plt.xticks(np.arange(0, 6, step=1))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend()
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_disp.eps", bbox_inches='tight', pad_inches=0.1)

Annotated version of examples/simple_shear_cyclic_loading/plot_reaction_force.py

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# import matplotlib.ticker as mtick
time1, forceX1, forceY1 = np.loadtxt('Reaction_force.hist',
delimiter='\t', unpack=True)
#time2, forceX2, forceY2 = np.loadtxt('Reaction_force_timestep_10-5.hist',
# delimiter='\t', unpack=True)
font = {'size': 18}
matplotlib.rc('font', **font)
labels = []
labels.append("Loading")
labels.append("Unloading")
labels.append("Rev. loading")
labels.append("Rev. unloading")
symbols = ['-k', '--*m', '-.b^', '--k', '-rD', 'bx']
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:16], forceX1[:16], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(time1[15:], forceX1[15:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Pseudo time')
plt.ylabel('Reaction force (kN)', multialignment='center')
plt.xlim(0,0.03)
# plt.xticks(np.arange(0, 6, step=1))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend(fontsize=16)
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_time.eps", bbox_inches='tight', pad_inches=0.1)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:16], forceX1[:16], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(30.0e-3 - time1[15:], forceX1[15:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Displacement (mm)')
plt.ylabel('Reaction force (kN)', multialignment='center')
plt.xlim(0.0,0.02)
#plt.xticks(np.arange(-0.02, 0.021, step=0.01))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend(fontsize=16)
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_disp.eps", bbox_inches='tight', pad_inches=0.1)

Annotated version of examples/simple_tension_cyclic_loading/plot_reaction_force.py

import matplotlib
import numpy as np
import matplotlib.pyplot as plt
# import matplotlib.ticker as mtick
time1, forceX1, forceY1 = np.loadtxt('Reaction_force.hist',
delimiter='\t', unpack=True)
#time2, forceX2, forceY2 = np.loadtxt('Reaction_force_timestep_10-5.hist',
# delimiter='\t', unpack=True)
font = {'size': 18}
matplotlib.rc('font', **font)
labels = []
labels.append("Loading")
labels.append("Unloading")
symbols = ['-k', '--*m', '-.b^', ':gx', '-rD', 'bx']
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:25], forceY1[:25], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(time1[25:], forceY1[25:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Pseudo time')
plt.ylabel('Reaction force (kN)', multialignment='center')
#plt.xlim(0,0.05)
# plt.xticks(np.arange(0, 6, step=1))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend()
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_time.eps", bbox_inches='tight', pad_inches=0.1)
fig = plt.figure()
ax = fig.add_subplot(111)
plt.plot(time1[:25], forceY1[:25], symbols[0],
linewidth=2.0, label=labels[0], fillstyle='none', markersize=8)
plt.plot(14.0e-3 - time1[25:], forceY1[25:], symbols[1],
linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
#plt.plot(time2, forceX2, symbols[1],
# linewidth=2.0, label=labels[1], fillstyle='none', markersize=8)
# plt.ticklabel_format(axis='y', style='sci', scilimits=(0, 3))
# ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
plt.xlabel('Displacement (mm)')
plt.ylabel('Reaction force (kN)', multialignment='center')
#plt.xlim(0,0.05)
# plt.xticks(np.arange(0, 6, step=1))
#plt.ylim(0.0, 1.0)
# plt.legend( bbox_to_anchor=(2.4, 1.1),prop={'size':10})
plt.grid()
plt.legend()
# ax.set_aspect('equal')
# plt.subplots_adjust(top = 0.95, bottom = 0.125,
# right = 0.95, left = 0.125, hspace = 0, wspace = 0)
plt.savefig("Reaction_force_history_disp.eps", bbox_inches='tight', pad_inches=0.1)

Annotated version of main.cc

  /* ---------------------------------------------------------------------
  *
  * Copyright (C) 2006 - 2020 by the deal.II authors
  *
  * 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.md at
  * the top level directory of deal.II.
  *
  * ---------------------------------------------------------------------
 
  *
  * Author: Tao Jin
  * University of Ottawa, Ottawa, Ontario, Canada
  * April. 2024
  *
  * How to cite:
  * Jin T, Li Z, Chen K. A novel phase-field monolithic scheme for brittle crack
  * propagation based on the limited-memory BFGS method with adaptive mesh refinement.
  * Int J Numer Methods Eng. 2024;e7572. doi: 10.1002/nme.7572
  */
 
  /* A monolithic scheme based on the L-BFGS method to solve the phase-field crack problem
  * 1. The phase-field formulation itself is based on "A phase field model for rate-independent
  * crack propagation - Robust algorithmic implementation based on operator splits"
  * by Christian Miehe , Martina Hofacker, Fabian Welschinger
  * 2. This code implements a monolithic approach. The phase-field irreversibility
  * is enforced through the history field Phi_0^+ and the viscosity parameter.
  * 3. Using TBB for stiffness assembly and Gauss point calculation.
  * 4. Using adaptive mesh refinement.
  */
 
  #include <deal.II/grid/tria.h>
  #include <deal.II/grid/grid_generator.h>
  #include <deal.II/grid/grid_refinement.h>
  #include <deal.II/grid/grid_out.h>
  #include <deal.II/grid/grid_in.h>
  #include <deal.II/grid/manifold_lib.h>
 
  #include <deal.II/dofs/dof_handler.h>
  #include <deal.II/dofs/dof_tools.h>
  #include <deal.II/dofs/dof_renumbering.h>
 
  #include <deal.II/fe/fe_values.h>
  #include <deal.II/fe/fe_system.h>
  #include <deal.II/fe/fe_q.h>
  #include <deal.II/fe/fe_dgp_monomial.h>
  #include <deal.II/fe/mapping_q_eulerian.h>
 
  #include <deal.II/base/timer.h>
  #include <deal.II/base/quadrature_point_data.h>
  #include <deal.II/base/parameter_handler.h>
 
  #include <deal.II/lac/affine_constraints.h>
  #include <deal.II/lac/vector.h>
  #include <deal.II/lac/full_matrix.h>
  #include <deal.II/lac/sparse_matrix.h>
  #include <deal.II/lac/dynamic_sparsity_pattern.h>
  #include <deal.II/lac/block_sparse_matrix.h>
  #include <deal.II/lac/block_vector.h>
 
 
  #include <deal.II/numerics/vector_tools.h>
  #include <deal.II/numerics/matrix_tools.h>
  #include <deal.II/numerics/data_out.h>
 
  #include <deal.II/lac/solver_cg.h>
  #include <deal.II/lac/precondition.h>
  #include <deal.II/lac/linear_operator.h>
  #include <deal.II/lac/packaged_operation.h>
  #include <deal.II/lac/precondition_selector.h>
  #include <deal.II/lac/solver_selector.h>
  #include <deal.II/lac/sparse_direct.h>
 
  #include <deal.II/numerics/error_estimator.h>
 
  #include <deal.II/physics/elasticity/standard_tensors.h>
 
  #include <deal.II/base/quadrature_point_data.h>
 
  #include <deal.II/grid/grid_tools.h>
 
  #include <deal.II/base/work_stream.h>
 
  #include <deal.II/numerics/solution_transfer.h>
 
  #include <fstream>
  #include <iostream>
  #include <deal.II/base/logstream.h>
 
  #include "SpectrumDecomposition.h"
  #include "Utilities.h"
 
  namespace PhaseField
  {
  using namespace dealii;
 
linear_operator
LinearOperator< Range, Domain, Payload > linear_operator(const OperatorExemplar &, const Matrix &)
Definition linear_operator.h:1490

body force

  template <int dim>
  void right_hand_side(const std::vector<Point<dim>> &points,
  std::vector<Tensor<1, dim>> & values,
  const double fx,
  const double fy,
  const double fz)
  {
  Assert(values.size() == points.size(),
  ExcDimensionMismatch(values.size(), points.size()));
  Assert(dim >= 2, ExcNotImplemented());
 
  for (unsigned int point_n = 0; point_n < points.size(); ++point_n)
  {
  if (dim == 2)
  {
  values[point_n][0] = fx;
  values[point_n][1] = fy;
  }
  else
  {
  values[point_n][0] = fx;
  values[point_n][1] = fy;
  values[point_n][2] = fz;
  }
  }
  }
 
  double degradation_function(const double d)
  {
  return (1.0 - d) * (1.0 - d);
  }
 
  double degradation_function_derivative(const double d)
  {
  return 2.0 * (d - 1.0);
  }
 
  double degradation_function_2nd_order_derivative(const double d)
  {
  (void) d;
  return 2.0;
  }
 
  namespace Parameters
  {
  struct Scenario
  {
  unsigned int m_scenario;
  std::string m_logfile_name;
  bool m_output_iteration_history;
  std::string m_type_nonlinear_solver;
  std::string m_type_line_search;
  std::string m_type_linear_solver;
  std::string m_refinement_strategy;
  unsigned int m_LBFGS_m;
  unsigned int m_global_refine_times;
  unsigned int m_local_prerefine_times;
  unsigned int m_max_adaptive_refine_times;
  int m_max_allowed_refinement_level;
  double m_phasefield_refine_threshold;
  double m_allowed_max_h_l_ratio;
  unsigned int m_total_material_regions;
  std::string m_material_file_name;
  int m_reaction_force_face_id;
 
  static void declare_parameters(ParameterHandler &prm);
  void parse_parameters(ParameterHandler &prm);
  };
 
  void Scenario::declare_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Scenario");
  {
  prm.declare_entry("Scenario number",
  "1",
  Patterns::Integer(0),
  "Geometry, loading and boundary conditions scenario");
 
  prm.declare_entry("Log file name",
  "Output.log",
  Patterns::FileName(Patterns::FileName::input),
  "Name of the file for log");
 
  prm.declare_entry("Output iteration history",
  "yes",
  Patterns::Selection("yes|no"),
  "Shall we write iteration history to the log file?");
 
  prm.declare_entry("Nonlinear solver type",
  "Newton",
  Patterns::Selection("Newton|BFGS|LBFGS"),
  "Type of solver used to solve the nonlinear system");
 
  prm.declare_entry("Line search type",
  "GradientBased",
  Patterns::Selection("GradientBased|StrongWolfe"),
  "Type of line search method, the gradient-based method "
  "should be preferred since it is generally faster");
 
  prm.declare_entry("Linear solver type",
  "Direct",
  Patterns::Selection("Direct|CG"),
  "Type of solver used to solve the linear system B0");
 
  prm.declare_entry("Mesh refinement strategy",
  "adaptive-refine",
  Patterns::Selection("pre-refine|adaptive-refine"),
  "Mesh refinement strategy: pre-refine or adaptive-refine");
 
  prm.declare_entry("LBFGS m",
  "40",
  Patterns::Integer(0),
  "Number of vectors used for LBFGS");
 
  prm.declare_entry("Global refinement times",
  "0",
  Patterns::Integer(0),
  "Global refinement times (across the entire domain)");
 
  prm.declare_entry("Local prerefinement times",
  "0",
  Patterns::Integer(0),
  "Local pre-refinement times (assume crack path is known a priori), "
  "only refine along the crack path.");
 
  prm.declare_entry("Max adaptive refinement times",
  "100",
  Patterns::Integer(0),
  "Maximum number of adaptive refinement times allowed in each step");
 
  prm.declare_entry("Max allowed refinement level",
  "100",
  Patterns::Integer(0),
  "Maximum allowed cell refinement level");
 
  prm.declare_entry("Phasefield refine threshold",
  "0.8",
  Patterns::Double(),
  "Phasefield-based refinement threshold value");
 
  prm.declare_entry("Allowed max hl ratio",
  "0.25",
  Patterns::Double(),
  "Allowed maximum ratio between mesh size h and length scale l");
 
  prm.declare_entry("Material regions",
  "1",
  Patterns::Integer(0),
  "Number of material regions");
 
  prm.declare_entry("Material data file",
  "1",
  Patterns::FileName(Patterns::FileName::input),
  "Material data file");
 
  prm.declare_entry("Reaction force face ID",
  "1",
  Patterns::Integer(),
  "Face id where reaction forces should be calculated "
  "(negative integer means not to calculate reaction force)");
  }
  prm.leave_subsection();
  }
 
  void Scenario::parse_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Scenario");
  {
  m_scenario = prm.get_integer("Scenario number");
  m_logfile_name = prm.get("Log file name");
  m_output_iteration_history = prm.get_bool("Output iteration history");
  m_type_nonlinear_solver = prm.get("Nonlinear solver type");
  m_type_line_search = prm.get("Line search type");
  m_type_linear_solver = prm.get("Linear solver type");
  m_refinement_strategy = prm.get("Mesh refinement strategy");
  m_LBFGS_m = prm.get_integer("LBFGS m");
  m_global_refine_times = prm.get_integer("Global refinement times");
  m_local_prerefine_times = prm.get_integer("Local prerefinement times");
  m_max_adaptive_refine_times = prm.get_integer("Max adaptive refinement times");
  m_max_allowed_refinement_level = prm.get_integer("Max allowed refinement level");
  m_phasefield_refine_threshold = prm.get_double("Phasefield refine threshold");
  m_allowed_max_h_l_ratio = prm.get_double("Allowed max hl ratio");
  m_total_material_regions = prm.get_integer("Material regions");
  m_material_file_name = prm.get("Material data file");
  m_reaction_force_face_id = prm.get_integer("Reaction force face ID");
  }
  prm.leave_subsection();
  }
 
  struct FESystem
  {
  unsigned int m_poly_degree;
  unsigned int m_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(0),
  "Phase field polynomial order");
 
  prm.declare_entry("Quadrature order",
  "2",
  Patterns::Integer(0),
  "Gauss quadrature order");
  }
  prm.leave_subsection();
  }
 
  void FESystem::parse_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Finite element system");
  {
  m_poly_degree = prm.get_integer("Polynomial degree");
  m_quad_order = prm.get_integer("Quadrature order");
  }
  prm.leave_subsection();
  }
 
ParameterHandler::enter_subsection
void enter_subsection(const std::string &subsection, const bool create_path_if_needed=true)
Definition parameter_handler.cc:1064
ParameterHandler::get_integer
long int get_integer(const std::string &entry_string) const
Definition parameter_handler.cc:1160
ParameterHandler::get_bool
bool get_bool(const std::string &entry_name) const
Definition parameter_handler.cc:1248
ParameterHandler::declare_entry
void declare_entry(const std::string &entry, const std::string &default_value, const Patterns::PatternBase &pattern=Patterns::Anything(), const std::string &documentation="", const bool has_to_be_set=false)
Definition parameter_handler.cc:889
ParameterHandler::get
std::string get(const std::string &entry_string) const
Definition parameter_handler.cc:1122
ParameterHandler::get_double
double get_double(const std::string &entry_name) const
Definition parameter_handler.cc:1203
Point
Definition point.h:111
Physics::Elasticity::Kinematics::d
SymmetricTensor< 2, dim, Number > d(const Tensor< 2, dim, Number > &F, const Tensor< 2, dim, Number > &dF_dt)

body force (N/m^3)

  struct BodyForce
  {
  double m_x_component;
  double m_y_component;
  double m_z_component;
 
  static void declare_parameters(ParameterHandler &prm);
 
  void parse_parameters(ParameterHandler &prm);
  };
 
  void BodyForce::declare_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Body force");
  {
  prm.declare_entry("Body force x component",
  "0.0",
  Patterns::Double(),
  "Body force x-component (N/m^3)");
 
  prm.declare_entry("Body force y component",
  "0.0",
  Patterns::Double(),
  "Body force y-component (N/m^3)");
 
  prm.declare_entry("Body force z component",
  "0.0",
  Patterns::Double(),
  "Body force z-component (N/m^3)");
  }
  prm.leave_subsection();
  }
 
  void BodyForce::parse_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Body force");
  {
  m_x_component = prm.get_double("Body force x component");
  m_y_component = prm.get_double("Body force y component");
  m_z_component = prm.get_double("Body force z component");
  }
  prm.leave_subsection();
  }
 
  struct NonlinearSolver
  {
  unsigned int m_max_iterations_NR;
  unsigned int m_max_iterations_BFGS;
  bool m_relative_residual;
 
  double m_tol_u_residual;
  double m_tol_d_residual;
  double m_tol_u_incr;
  double m_tol_d_incr;
 
  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("Max iterations BFGS",
  "20",
  Patterns::Integer(0),
  "Number of BFGS iterations allowed");
 
  prm.declare_entry("Relative residual",
  "yes",
  Patterns::Selection("yes|no"),
  "Shall we use relative residual for convergence?");
 
  prm.declare_entry("Tolerance displacement residual",
  "1.0e-9",
  Patterns::Double(0.0),
  "Displacement residual tolerance");
 
  prm.declare_entry("Tolerance phasefield residual",
  "1.0e-9",
  Patterns::Double(0.0),
  "Phasefield residual tolerance");
 
  prm.declare_entry("Tolerance displacement increment",
  "1.0e-9",
  Patterns::Double(0.0),
  "Displacement increment tolerance");
 
  prm.declare_entry("Tolerance phasefield increment",
  "1.0e-9",
  Patterns::Double(0.0),
  "Phasefield increment tolerance");
  }
  prm.leave_subsection();
  }
 
  void NonlinearSolver::parse_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Nonlinear solver");
  {
  m_max_iterations_NR = prm.get_integer("Max iterations Newton-Raphson");
  m_max_iterations_BFGS = prm.get_integer("Max iterations BFGS");
  m_relative_residual = prm.get_bool("Relative residual");
 
  m_tol_u_residual = prm.get_double("Tolerance displacement residual");
  m_tol_d_residual = prm.get_double("Tolerance phasefield residual");
  m_tol_u_incr = prm.get_double("Tolerance displacement increment");
  m_tol_d_incr = prm.get_double("Tolerance phasefield increment");
  }
  prm.leave_subsection();
  }
 
  struct TimeInfo
  {
  double m_end_time;
  std::string m_time_file_name;
 
  static void declare_parameters(ParameterHandler &prm);
 
  void parse_parameters(ParameterHandler &prm);
  };
 
  void TimeInfo::declare_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Time");
  {
  prm.declare_entry("End time", "1", Patterns::Double(), "End time");
 
  prm.declare_entry("Time data file",
  "1",
  Patterns::FileName(Patterns::FileName::input),
  "Time data file");
  }
  prm.leave_subsection();
  }
 
  void TimeInfo::parse_parameters(ParameterHandler &prm)
  {
  prm.enter_subsection("Time");
  {
  m_end_time = prm.get_double("End time");
  m_time_file_name = prm.get("Time data file");
  }
  prm.leave_subsection();
  }
 
  struct AllParameters : public Scenario,
  public FESystem,
  public BodyForce,
  public NonlinearSolver,
  public TimeInfo
  {
  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)
  {
  Scenario::declare_parameters(prm);
  FESystem::declare_parameters(prm);
  BodyForce::declare_parameters(prm);
  NonlinearSolver::declare_parameters(prm);
  TimeInfo::declare_parameters(prm);
  }
 
  void AllParameters::parse_parameters(ParameterHandler &prm)
  {
  Scenario::parse_parameters(prm);
  FESystem::parse_parameters(prm);
  BodyForce::parse_parameters(prm);
  NonlinearSolver::parse_parameters(prm);
  TimeInfo::parse_parameters(prm);
  }
  } // namespace Parameters
 
  class Time
  {
  public:
  Time(const double time_end)
  : m_timestep(0)
  , m_time_current(0.0)
  , m_time_end(time_end)
  , m_delta_t(0.0)
  , m_magnitude(1.0)
  {}
 
  virtual ~Time() = default;
 
  double current() const
  {
  return m_time_current;
  }
  double end() const
  {
  return m_time_end;
  }
  double get_delta_t() const
  {
  return m_delta_t;
  }
  double get_magnitude() const
  {
  return m_magnitude;
  }
  unsigned int get_timestep() const
  {
  return m_timestep;
  }
  void increment(std::vector<std::array<double, 4>> time_table)
  {
  double t_1, t_delta, t_magnitude;
  for (auto & time_group : time_table)
  {
  t_1 = time_group[1];
  t_delta = time_group[2];
  t_magnitude = time_group[3];
 
  if (m_time_current < t_1 - 1.0e-6*t_delta)
  {
  m_delta_t = t_delta;
  m_magnitude = t_magnitude;
  break;
  }
  }
 
  m_time_current += m_delta_t;
  ++m_timestep;
  }
 
  private:
  unsigned int m_timestep;
  double m_time_current;
  const double m_time_end;
  double m_delta_t;
  double m_magnitude;
  };
 
  template <int dim>
  class LinearIsotropicElasticityAdditiveSplit
  {
  public:
  LinearIsotropicElasticityAdditiveSplit(const double lame_lambda,
  const double lame_mu,
  const double residual_k,
  const double length_scale,
  const double viscosity,
  const double gc)
  : m_lame_lambda(lame_lambda)
  , m_lame_mu(lame_mu)
  , m_residual_k(residual_k)
  , m_length_scale(length_scale)
  , m_eta(viscosity)
  , m_gc(gc)
  , m_phase_field_value(0.0)
  , m_grad_phasefield(Tensor<1, dim>())
  , m_strain(SymmetricTensor<2, dim>())
  , m_stress(SymmetricTensor<2, dim>())
  , m_stress_positive(SymmetricTensor<2, dim>())
  , m_mechanical_C(SymmetricTensor<4, dim>())
  , m_strain_energy_positive(0.0)
  , m_strain_energy_negative(0.0)
  , m_strain_energy_total(0.0)
  , m_crack_energy_dissipation(0.0)
  {
  Assert( ( lame_lambda / (2*(lame_lambda + lame_mu)) <= 0.5)
  & ( lame_lambda / (2*(lame_lambda + lame_mu)) >=-1.0),
  ExcInternalError() );
  }
 
  const SymmetricTensor<4, dim> & get_mechanical_C() const
  {
  return m_mechanical_C;
  }
 
  const SymmetricTensor<2, dim> & get_cauchy_stress() const
  {
  return m_stress;
  }
 
  const SymmetricTensor<2, dim> & get_cauchy_stress_positive() const
  {
  return m_stress_positive;
  }
 
  double get_positive_strain_energy() const
  {
  return m_strain_energy_positive;
  }
 
  double get_negative_strain_energy() const
  {
  return m_strain_energy_negative;
  }
 
  double get_total_strain_energy() const
  {
  return m_strain_energy_total;
  }
 
  double get_crack_energy_dissipation() const
  {
  return m_crack_energy_dissipation;
  }
 
  double get_phase_field_value() const
  {
  return m_phase_field_value;
  }
 
  const Tensor<1, dim> get_phase_field_gradient() const
  {
  return m_grad_phasefield;
  }
 
  void update_material_data(const SymmetricTensor<2, dim> & strain,
  const double phase_field_value,
  const Tensor<1, dim> & grad_phasefield,
  const double phase_field_value_previous_step,
  const double delta_time)
  {
  m_strain = strain;
  m_phase_field_value = phase_field_value;
  m_grad_phasefield = grad_phasefield;
  Vector<double> eigenvalues(dim);
  std::vector<Tensor<1, dim>> eigenvectors(dim);
  usr_spectrum_decomposition::spectrum_decomposition<dim>(m_strain,
  eigenvalues,
  eigenvectors);
 
  SymmetricTensor<2, dim> strain_positive, strain_negative;
  strain_positive = usr_spectrum_decomposition::positive_tensor(eigenvalues, eigenvectors);
  strain_negative = usr_spectrum_decomposition::negative_tensor(eigenvalues, eigenvectors);
 
  SymmetricTensor<4, dim> projector_positive, projector_negative;
  usr_spectrum_decomposition::positive_negative_projectors(eigenvalues,
  eigenvectors,
  projector_positive,
  projector_negative);
 
  SymmetricTensor<2, dim> stress_positive, stress_negative;
  const double degradation = degradation_function(m_phase_field_value) + m_residual_k;
  const double I_1 = trace(m_strain);
  stress_positive = m_lame_lambda * usr_spectrum_decomposition::positive_ramp_function(I_1)
  * Physics::Elasticity::StandardTensors<dim>::I
  + 2 * m_lame_mu * strain_positive;
  stress_negative = m_lame_lambda * usr_spectrum_decomposition::negative_ramp_function(I_1)
  * Physics::Elasticity::StandardTensors<dim>::I
  + 2 * m_lame_mu * strain_negative;
 
  m_stress = degradation * stress_positive + stress_negative;
  m_stress_positive = stress_positive;
 
  SymmetricTensor<4, dim> C_positive, C_negative;
  C_positive = m_lame_lambda * usr_spectrum_decomposition::heaviside_function(I_1)
  * Physics::Elasticity::StandardTensors<dim>::IxI
  + 2 * m_lame_mu * projector_positive;
  C_negative = m_lame_lambda * usr_spectrum_decomposition::heaviside_function(-I_1)
  * Physics::Elasticity::StandardTensors<dim>::IxI
  + 2 * m_lame_mu * projector_negative;
  m_mechanical_C = degradation * C_positive + C_negative;
 
  m_strain_energy_positive = 0.5 * m_lame_lambda * usr_spectrum_decomposition::positive_ramp_function(I_1)
  * usr_spectrum_decomposition::positive_ramp_function(I_1)
  + m_lame_mu * strain_positive * strain_positive;
 
  m_strain_energy_negative = 0.5 * m_lame_lambda * usr_spectrum_decomposition::negative_ramp_function(I_1)
  * usr_spectrum_decomposition::negative_ramp_function(I_1)
  + m_lame_mu * strain_negative * strain_negative;
 
  m_strain_energy_total = degradation * m_strain_energy_positive + m_strain_energy_negative;
 
  m_crack_energy_dissipation = m_gc * ( 0.5 / m_length_scale * m_phase_field_value * m_phase_field_value
  + 0.5 * m_length_scale * m_grad_phasefield * m_grad_phasefield)
end
*  iterator end()
Definition ALE_Finite_Strain_Plasticity.h:10398
trace
constexpr Number trace(const SymmetricTensor< 2, dim2, Number > &)

the term due to viscosity regularization

  + (m_phase_field_value - phase_field_value_previous_step)
  * (m_phase_field_value - phase_field_value_previous_step)
  * 0.5 * m_eta / delta_time;

(void)delta_time; (void)phase_field_value_previous_step;

  }
 
  private:
  const double m_lame_lambda;
  const double m_lame_mu;
  const double m_residual_k;
  const double m_length_scale;
  const double m_eta;
  const double m_gc;
  double m_phase_field_value;
  Tensor<1, dim> m_grad_phasefield;
  SymmetricTensor<2, dim> m_strain;
  SymmetricTensor<2, dim> m_stress;
  SymmetricTensor<2, dim> m_stress_positive;
  SymmetricTensor<4, dim> m_mechanical_C;
  double m_strain_energy_positive;
  double m_strain_energy_negative;
  double m_strain_energy_total;
  double m_crack_energy_dissipation;
  };
 
 
  template <int dim>
  class PointHistory
  {
  public:
  PointHistory()
  : m_length_scale(0.0)
  , m_gc(0.0)
  , m_viscosity(0.0)
  , m_history_max_positive_strain_energy(0.0)
  {}
 
  virtual ~PointHistory() = default;
 
  void setup_lqp(const double lame_lambda,
  const double lame_mu,
  const double length_scale,
  const double gc,
  const double viscosity,
  const double residual_k)
  {
  m_material =
  std::make_shared<LinearIsotropicElasticityAdditiveSplit<dim>>(lame_lambda,
  lame_mu,
  residual_k,
  length_scale,
  viscosity,
  gc);
  m_history_max_positive_strain_energy = 0.0;
  m_length_scale = length_scale;
  m_gc = gc;
  m_viscosity = viscosity;
 
  update_field_values(SymmetricTensor<2, dim>(), 0.0, Tensor<1, dim>(), 0.0, 1.0);
  }
 
  void update_field_values(const SymmetricTensor<2, dim> & strain,
  const double phase_field_value,
  const Tensor<1, dim> & grad_phasefield,
  const double phase_field_value_previous_step,
  const double delta_time)
  {
  m_material->update_material_data(strain, phase_field_value, grad_phasefield,
  phase_field_value_previous_step, delta_time);
  }
 
  void update_history_variable()
  {
  double current_positive_strain_energy = m_material->get_positive_strain_energy();
  m_history_max_positive_strain_energy = std::fmax(m_history_max_positive_strain_energy,
  current_positive_strain_energy);
  }
 
  double get_current_positive_strain_energy() const
  {
  return m_material->get_positive_strain_energy();
  }
 
  const SymmetricTensor<4, dim> & get_mechanical_C() const
  {
  return m_material->get_mechanical_C();
  }
 
  const SymmetricTensor<2, dim> & get_cauchy_stress() const
  {
  return m_material->get_cauchy_stress();
  }
 
  const SymmetricTensor<2, dim> & get_cauchy_stress_positive() const
  {
  return m_material->get_cauchy_stress_positive();
  }
 
  double get_total_strain_energy() const
  {
  return m_material->get_total_strain_energy();
  }
 
  double get_crack_energy_dissipation() const
  {
  return m_material->get_crack_energy_dissipation();
  }
 
  double get_phase_field_value() const
  {
  return m_material->get_phase_field_value();
  }
 
  const Tensor<1, dim> get_phase_field_gradient() const
  {
  return m_material->get_phase_field_gradient();
  }
 
  double get_history_max_positive_strain_energy() const
  {
  return m_history_max_positive_strain_energy;
  }
 
  double get_length_scale() const
  {
  return m_length_scale;
  }
 
  double get_critical_energy_release_rate() const
  {
  return m_gc;
  }
 
  double get_viscosity() const
  {
  return m_viscosity;
  }
  private:
  std::shared_ptr<LinearIsotropicElasticityAdditiveSplit<dim>> m_material;
  double m_length_scale;
  double m_gc;
  double m_viscosity;
  double m_history_max_positive_strain_energy;
  };
 
  template <int dim>
  class PhaseFieldMonolithicSolve
  {
  public:
  PhaseFieldMonolithicSolve(const std::string &input_file);
 
  virtual ~PhaseFieldMonolithicSolve() = default;
  void run();
 
  private:
  struct PerTaskData_ASM;
  struct ScratchData_ASM;
 
  struct PerTaskData_ASM_RHS_BFGS;
  struct ScratchData_ASM_RHS_BFGS;
 
  struct PerTaskData_UQPH;
  struct ScratchData_UQPH;
 
  Parameters::AllParameters m_parameters;
  Triangulation<dim> m_triangulation;
 
  CellDataStorage<typename Triangulation<dim>::cell_iterator,
  PointHistory<dim>>
  m_quadrature_point_history;
 
  Time m_time;
  std::ofstream m_logfile;
  mutable TimerOutput m_timer;
 
  DoFHandler<dim> m_dof_handler;
  FESystem<dim> m_fe;
  const unsigned int m_dofs_per_cell;
  const FEValuesExtractors::Vector m_u_fe;
  const FEValuesExtractors::Scalar m_d_fe;
 
  static const unsigned int m_n_blocks = 2;
  static const unsigned int m_n_components = dim + 1;
  static const unsigned int m_first_u_component = 0;
  static const unsigned int m_d_component = dim;
 
  enum
  {
  m_u_dof = 0,
  m_d_dof = 1
  };
 
  std::vector<types::global_dof_index> m_dofs_per_block;
 
  const QGauss<dim> m_qf_cell;
  const QGauss<dim - 1> m_qf_face;
  const unsigned int m_n_q_points;
 
  double m_vol_reference;
 
  AffineConstraints<double> m_constraints;
  BlockSparsityPattern m_sparsity_pattern;
  BlockSparseMatrix<double> m_tangent_matrix;
  BlockVector<double> m_system_rhs;
  BlockVector<double> m_solution;
  SparseDirectUMFPACK m_A_direct;
 
 
  std::map<unsigned int, std::vector<double>> m_material_data;
 
  std::vector<std::pair<double, std::vector<double>>> m_history_reaction_force;
  std::vector<std::pair<double, std::array<double, 3>>> m_history_energy;
 
 
  struct Errors
  {
  Errors()
  : m_norm(1.0)
  , m_u(1.0)
  , m_d(1.0)
  {}
 
  void reset()
  {
  m_norm = 1.0;
  m_u = 1.0;
  m_d = 1.0;
  }
 
  void normalize(const Errors &rhs)
  {
  if (rhs.m_norm != 0.0)
  m_norm /= rhs.m_norm;
  if (rhs.m_u != 0.0)
  m_u /= rhs.m_u;
  if (rhs.m_d != 0.0)
  m_d /= rhs.m_d;
  }
 
  double m_norm, m_u, m_d;
  };
 
  Errors m_error_residual, m_error_residual_0, m_error_residual_norm, m_error_update,
  m_error_update_0, m_error_update_norm;
 
  void get_error_residual(Errors &error_residual);
  void get_error_update(const BlockVector<double> &newton_update,
  Errors & error_update);
 
  void make_grid();
  void make_grid_case_1();
  void make_grid_case_2();
  void make_grid_case_3();
  void make_grid_case_4();
  void make_grid_case_5();
  void make_grid_case_6();
  void make_grid_case_7();
  void make_grid_case_8();
  void make_grid_case_9();
  void make_grid_case_11();
 
  void setup_system();
 
  void determine_component_extractors();
 
  void make_constraints(const unsigned int it_nr);
 
  void assemble_system_newton(const BlockVector<double> & solution_old);
 
  void assemble_system_B0(const BlockVector<double> & solution_old);
 
  void assemble_system_newton_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data) const;
 
  void assemble_system_B0_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data) const;
 
  void assemble_system_rhs_BFGS_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM_RHS_BFGS & scratch,
  PerTaskData_ASM_RHS_BFGS & data) const;
 
  void assemble_system_rhs_BFGS(const BlockVector<double> & solution_old,
  BlockVector<double> & system_rhs);
 
  void assemble_system_rhs_BFGS_parallel(const BlockVector<double> & solution_old,
  BlockVector<double> & system_rhs);
 
  bool solve_nonlinear_timestep_newton(BlockVector<double> &solution_delta);
 
  void solve_nonlinear_timestep_BFGS(BlockVector<double> &solution_delta);
 
  void solve_nonlinear_timestep_LBFGS(BlockVector<double> &solution_delta,
  BlockVector<double> & LBFGS_update_refine);
 
  double line_search_stepsize_strong_wolfe(const double phi_0,
  const double phi_0_prime,
  const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta);
 
  double line_search_stepsize_gradient_based(const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta);
 
  double line_search_zoom_strong_wolfe(double phi_low, double phi_low_prime, double alpha_low,
  double phi_high, double phi_high_prime, double alpha_high,
  double phi_0, double phi_0_prime, const BlockVector<double> & BFGS_p_vector,
  double c1, double c2, unsigned int max_iter,
  const BlockVector<double> & solution_delta);
 
  double line_search_interpolation_cubic(const double alpha_0, const double phi_0, const double phi_0_prime,
  const double alpha_1, const double phi_1, const double phi_1_prime);
 
  std::pair<double, double> calculate_phi_and_phi_prime(const double alpha,
  const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta);
 
  std::vector<double> solve_linear_system(BlockVector<double> &newton_update);
 
  void LBFGS_B0(BlockVector<double> & LBFGS_r_vector,
  BlockVector<double> & LBFGS_q_vector);
 
  void update_history_field_step();
 
  void output_results() const;
 
  void setup_qph();
 
  void update_qph_incremental(const BlockVector<double> &solution_delta,
  const BlockVector<double> &solution_old,
  const bool is_print);
 
  void update_qph_incremental_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_UQPH & scratch,
  PerTaskData_UQPH & data);
 
  void copy_local_to_global_UQPH(const PerTaskData_UQPH & /*data*/)
  {}
 
  BlockVector<double>
  get_total_solution(const BlockVector<double> &solution_delta) const;
 
DoFHandler::active_cell_iterator
typename ActiveSelector::active_cell_iterator active_cell_iterator
Definition dof_handler.h:447
data
std::vector< index_type > data
Definition mpi.cc:734
WorkStream::internal::tbb_no_coloring::run
void run(const Iterator &begin, const std_cxx20::type_identity_t< Iterator > &end, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length, const unsigned int chunk_size)
Definition work_stream.h:517

Should not make this function const

  void read_material_data(const std::string &data_file,
  const unsigned int total_material_regions);
 
  void read_time_data(const std::string &data_file,
  std::vector<std::array<double, 4>> & time_table);
 
  void print_conv_header_newton();
 
  void print_conv_header_BFGS();
 
  void print_conv_header_LBFGS();
 
  void print_parameter_information();
 
  void calculate_reaction_force(unsigned int face_ID);
 
  void write_history_data();
 
  double calculate_energy_functional() const;
 
  std::pair<double, double> calculate_total_strain_energy_and_crack_energy_dissipation() const;
 
  bool local_refine_and_solution_transfer(BlockVector<double> & solution_delta,
  BlockVector<double> & LBFGS_update_refine);
  }; // class PhaseFieldSplitSolve
 
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::get_error_residual(Errors &error_residual)
  {
  BlockVector<double> error_res(m_dofs_per_block);
 
  for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
  if (!m_constraints.is_constrained(i))
  error_res(i) = m_system_rhs(i);
 
  error_residual.m_norm = error_res.l2_norm();
  error_residual.m_u = error_res.block(m_u_dof).l2_norm();
  error_residual.m_d = error_res.block(m_d_dof).l2_norm();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::get_error_update(const BlockVector<double> &newton_update,
  Errors & error_update)
  {
  BlockVector<double> error_ud(m_dofs_per_block);
  for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i)
  if (!m_constraints.is_constrained(i))
  error_ud(i) = newton_update(i);
 
  error_update.m_norm = error_ud.l2_norm();
  error_update.m_u = error_ud.block(m_u_dof).l2_norm();
  error_update.m_d = error_ud.block(m_d_dof).l2_norm();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::read_material_data(const std::string &data_file,
  const unsigned int total_material_regions)
  {
  std::ifstream myfile (data_file);
 
  double lame_lambda, lame_mu, length_scale, gc, viscosity, residual_k;
  int material_region;
  double poisson_ratio;
  if (myfile.is_open())
  {
  m_logfile << "Reading material data file ..." << std::endl;
 
  while ( myfile >> material_region
  >> lame_lambda
  >> lame_mu
  >> length_scale
  >> gc
  >> viscosity
  >> residual_k)
  {
  m_material_data[material_region] = {lame_lambda,
  lame_mu,
  length_scale,
  gc,
  viscosity,
  residual_k};
  poisson_ratio = lame_lambda / (2*(lame_lambda + lame_mu));
  Assert( (poisson_ratio <= 0.5)&(poisson_ratio >=-1.0) , ExcInternalError());
 
  m_logfile << "\tRegion " << material_region << " : " << std::endl;
  m_logfile << "\t\tLame lambda = " << lame_lambda << std::endl;
  m_logfile << "\t\tLame mu = " << lame_mu << std::endl;
  m_logfile << "\t\tPoisson ratio = " << poisson_ratio << std::endl;
  m_logfile << "\t\tPhase field length scale (l) = " << length_scale << std::endl;
  m_logfile << "\t\tCritical energy release rate (gc) = " << gc << std::endl;
  m_logfile << "\t\tViscosity for regularization (eta) = " << viscosity << std::endl;
  m_logfile << "\t\tResidual_k (k) = " << residual_k << std::endl;
  }
 
  if (m_material_data.size() != total_material_regions)
  {
  m_logfile << "Material data file has " << m_material_data.size() << " rows. However, "
  << "the mesh has " << total_material_regions << " material regions."
  << std::endl;
  Assert(m_material_data.size() == total_material_regions,
  ExcDimensionMismatch(m_material_data.size(), total_material_regions));
  }
  myfile.close();
  }
  else
  {
  m_logfile << "Material data file : " << data_file << " not exist!" << std::endl;
  Assert(false, ExcMessage("Failed to read material data file"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::read_time_data(const std::string &data_file,
  std::vector<std::array<double, 4>> & time_table)
  {
  std::ifstream myfile (data_file);
 
  double t_0, t_1, delta_t, t_magnitude;
 
  if (myfile.is_open())
  {
  m_logfile << "Reading time data file ..." << std::endl;
 
  while ( myfile >> t_0
  >> t_1
  >> delta_t
  >> t_magnitude)
  {
  Assert( t_0 < t_1,
  ExcMessage("For each time pair, "
  "the start time should be smaller than the end time"));
  time_table.push_back({{t_0, t_1, delta_t, t_magnitude}});
  }
 
  Assert(std::fabs(t_1 - m_parameters.m_end_time) < 1.0e-9,
  ExcMessage("End time in time table is inconsistent with input data in parameters.prm"));
 
  Assert(time_table.size() > 0,
  ExcMessage("Time data file is empty."));
  myfile.close();
  }
  else
  {
  m_logfile << "Time data file : " << data_file << " not exist!" << std::endl;
  Assert(false, ExcMessage("Failed to read time data file"));
  }
 
  for (auto & time_group : time_table)
  {
  m_logfile << "\t\t"
  << time_group[0] << ",\t"
  << time_group[1] << ",\t"
  << time_group[2] << ",\t"
  << time_group[3] << std::endl;
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::setup_qph()
  {
  m_logfile << "\t\tSetting up quadrature point data ("
  << m_n_q_points
  << " points per cell)" << std::endl;
 
  m_quadrature_point_history.clear();
  for (auto const & cell : m_triangulation.active_cell_iterators())
  {
  m_quadrature_point_history.initialize(cell, m_n_q_points);
  }
 
  unsigned int material_id;
  double lame_lambda = 0.0;
  double lame_mu = 0.0;
  double length_scale = 0.0;
  double gc = 0.0;
  double viscosity = 0.0;
  double residual_k = 0.0;
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  material_id = cell->material_id();
  if (m_material_data.find(material_id) != m_material_data.end())
  {
  lame_lambda = m_material_data[material_id][0];
  lame_mu = m_material_data[material_id][1];
  length_scale = m_material_data[material_id][2];
  gc = m_material_data[material_id][3];
  viscosity = m_material_data[material_id][4];
  residual_k = m_material_data[material_id][5];
  }
  else
  {
  m_logfile << "Could not find material data for material id: " << material_id << std::endl;
  AssertThrow(false, ExcMessage("Could not find material data for material id."));
  }
 
  const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
  lqph[q_point]->setup_lqp(lame_lambda, lame_mu, length_scale,
  gc, viscosity, residual_k);
  }
  }
 
  template <int dim>
  BlockVector<double> PhaseFieldMonolithicSolve<dim>::get_total_solution(
  const BlockVector<double> &solution_delta) const
  {
  BlockVector<double> solution_total(m_solution);
  solution_total += solution_delta;
  return solution_total;
  }
 
  template <int dim>
  void
  PhaseFieldMonolithicSolve<dim>::update_qph_incremental(const BlockVector<double> &solution_delta,
  const BlockVector<double> &solution_old,
  const bool is_print)
  {
  m_timer.enter_subsection("Update QPH data");
  if (is_print && m_parameters.m_output_iteration_history)
  m_logfile << " UQPH " << std::flush;
 
  const BlockVector<double> solution_total(get_total_solution(solution_delta));
 
  const UpdateFlags uf_UQPH(update_values | update_gradients);
  PerTaskData_UQPH per_task_data_UQPH;
  ScratchData_UQPH scratch_data_UQPH(m_fe,
  m_qf_cell,
  uf_UQPH,
  solution_total,
  solution_old,
  m_time.get_delta_t());
 
  auto worker = [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_UQPH & scratch,
  PerTaskData_UQPH & data)
  {
  this->update_qph_incremental_one_cell(cell, scratch, data);
  };
 
  auto copier = [this](const PerTaskData_UQPH &data)
  {
  this->copy_local_to_global_UQPH(data);
  };
 
  WorkStream::run(
  m_dof_handler.begin_active(),
  m_dof_handler.end(),
  worker,
  copier,
  scratch_data_UQPH,
  per_task_data_UQPH);
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::PerTaskData_UQPH
  {
  void reset()
  {}
  };
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::ScratchData_UQPH
  {
  const BlockVector<double> & m_solution_UQPH;
 
  std::vector<SymmetricTensor<2, dim>> m_solution_symm_grads_u_cell;
  std::vector<double> m_solution_values_phasefield_cell;
  std::vector<Tensor<1, dim>> m_solution_grad_phasefield_cell;
 
  FEValues<dim> m_fe_values;
 
  const BlockVector<double>& m_solution_previous_step;
  std::vector<double> m_phasefield_previous_step_cell;
 
  const double m_delta_time;
 
  ScratchData_UQPH(const FiniteElement<dim> & fe_cell,
  const QGauss<dim> & qf_cell,
  const UpdateFlags uf_cell,
  const BlockVector<double> &solution_total,
  const BlockVector<double> &solution_old,
  const double delta_time)
  : m_solution_UQPH(solution_total)
  , m_solution_symm_grads_u_cell(qf_cell.size())
  , m_solution_values_phasefield_cell(qf_cell.size())
  , m_solution_grad_phasefield_cell(qf_cell.size())
  , m_fe_values(fe_cell, qf_cell, uf_cell)
  , m_solution_previous_step(solution_old)
  , m_phasefield_previous_step_cell(qf_cell.size())
  , m_delta_time(delta_time)
  {}
 
  ScratchData_UQPH(const ScratchData_UQPH &rhs)
  : m_solution_UQPH(rhs.m_solution_UQPH)
  , m_solution_symm_grads_u_cell(rhs.m_solution_symm_grads_u_cell)
  , m_solution_values_phasefield_cell(rhs.m_solution_values_phasefield_cell)
  , m_solution_grad_phasefield_cell(rhs.m_solution_grad_phasefield_cell)
  , m_fe_values(rhs.m_fe_values.get_fe(),
  rhs.m_fe_values.get_quadrature(),
  rhs.m_fe_values.get_update_flags())
  , m_solution_previous_step(rhs.m_solution_previous_step)
  , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
  , m_delta_time(rhs.m_delta_time)
  {}
 
  void reset()
  {
  const unsigned int n_q_points = m_solution_symm_grads_u_cell.size();
  for (unsigned int q = 0; q < n_q_points; ++q)
  {
  m_solution_symm_grads_u_cell[q] = 0.0;
  m_solution_values_phasefield_cell[q] = 0.0;
  m_solution_grad_phasefield_cell[q] = 0.0;
  m_phasefield_previous_step_cell[q] = 0.0;
  }
  }
  };
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::update_qph_incremental_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_UQPH & scratch,
  PerTaskData_UQPH & /*data*/)
  {
  scratch.reset();
 
  scratch.m_fe_values.reinit(cell);
 
  const std::vector<std::shared_ptr<PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  const FEValuesExtractors::Vector displacement(0);
 
  scratch.m_fe_values[m_u_fe].get_function_symmetric_gradients(
  scratch.m_solution_UQPH, scratch.m_solution_symm_grads_u_cell);
  scratch.m_fe_values[m_d_fe].get_function_values(
  scratch.m_solution_UQPH, scratch.m_solution_values_phasefield_cell);
  scratch.m_fe_values[m_d_fe].get_function_gradients(
  scratch.m_solution_UQPH, scratch.m_solution_grad_phasefield_cell);
 
  scratch.m_fe_values[m_d_fe].get_function_values(
  scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 
  for (const unsigned int q_point :
  scratch.m_fe_values.quadrature_point_indices())
  lqph[q_point]->update_field_values(scratch.m_solution_symm_grads_u_cell[q_point],
  scratch.m_solution_values_phasefield_cell[q_point],
  scratch.m_solution_grad_phasefield_cell[q_point],
  scratch.m_phasefield_previous_step_cell[q_point],
  scratch.m_delta_time);
  }
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::PerTaskData_ASM
  {
  FullMatrix<double> m_cell_matrix;
  Vector<double> m_cell_rhs;
  std::vector<types::global_dof_index> m_local_dof_indices;
 
  PerTaskData_ASM(const unsigned int dofs_per_cell)
  : m_cell_matrix(dofs_per_cell, dofs_per_cell)
  , m_cell_rhs(dofs_per_cell)
  , m_local_dof_indices(dofs_per_cell)
  {}
 
  void reset()
  {
  m_cell_matrix = 0.0;
  m_cell_rhs = 0.0;
  }
  };
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::PerTaskData_ASM_RHS_BFGS
  {
  Vector<double> m_cell_rhs;
  std::vector<types::global_dof_index> m_local_dof_indices;
 
  PerTaskData_ASM_RHS_BFGS(const unsigned int dofs_per_cell)
  : m_cell_rhs(dofs_per_cell)
  , m_local_dof_indices(dofs_per_cell)
  {}
 
  void reset()
  {
  m_cell_rhs = 0.0;
  }
  };
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::ScratchData_ASM
  {
  FEValues<dim> m_fe_values;
  FEFaceValues<dim> m_fe_face_values;
 
  std::vector<std::vector<double>> m_Nx_phasefield; // shape function values for phase-field
  std::vector<std::vector<Tensor<1, dim>>> m_grad_Nx_phasefield; // gradient of shape function values for phase field
 
  std::vector<std::vector<Tensor<1, dim>>> m_Nx_disp; // shape function values for displacement
  std::vector<std::vector<Tensor<2, dim>>> m_grad_Nx_disp; // gradient of shape function values for displacement
  std::vector<std::vector<SymmetricTensor<2, dim>>> m_symm_grad_Nx_disp; // symmetric gradient of shape function values for displacement
 
  const BlockVector<double>& m_solution_previous_step;
  std::vector<double> m_phasefield_previous_step_cell;
 
  ScratchData_ASM(const FiniteElement<dim> & fe_cell,
  const QGauss<dim> & qf_cell,
  const UpdateFlags uf_cell,
  const QGauss<dim - 1> & qf_face,
  const UpdateFlags uf_face,
  const BlockVector<double>& solution_old)
  : m_fe_values(fe_cell, qf_cell, uf_cell)
  , m_fe_face_values(fe_cell, qf_face, uf_face)
  , m_Nx_phasefield(qf_cell.size(),
  std::vector<double>(fe_cell.n_dofs_per_cell()))
  , m_grad_Nx_phasefield(qf_cell.size(),
  std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
  , m_Nx_disp(qf_cell.size(),
  std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
  , m_grad_Nx_disp(qf_cell.size(),
  std::vector<Tensor<2, dim>>(fe_cell.n_dofs_per_cell()))
  , m_symm_grad_Nx_disp(qf_cell.size(),
  std::vector<SymmetricTensor<2, dim>>(fe_cell.n_dofs_per_cell()))
  , m_solution_previous_step(solution_old)
  , m_phasefield_previous_step_cell(qf_cell.size())
  {}
 
  ScratchData_ASM(const ScratchData_ASM &rhs)
  : m_fe_values(rhs.m_fe_values.get_fe(),
  rhs.m_fe_values.get_quadrature(),
  rhs.m_fe_values.get_update_flags())
  , m_fe_face_values(rhs.m_fe_face_values.get_fe(),
  rhs.m_fe_face_values.get_quadrature(),
  rhs.m_fe_face_values.get_update_flags())
  , m_Nx_phasefield(rhs.m_Nx_phasefield)
  , m_grad_Nx_phasefield(rhs.m_grad_Nx_phasefield)
  , m_Nx_disp(rhs.m_Nx_disp)
  , m_grad_Nx_disp(rhs.m_grad_Nx_disp)
  , m_symm_grad_Nx_disp(rhs.m_symm_grad_Nx_disp)
  , m_solution_previous_step(rhs.m_solution_previous_step)
  , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
  {}
 
  void reset()
  {
  const unsigned int n_q_points = m_Nx_phasefield.size();
  const unsigned int n_dofs_per_cell = m_Nx_phasefield[0].size();
  for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
  {
  Assert(m_Nx_phasefield[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_grad_Nx_phasefield[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_symm_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  m_phasefield_previous_step_cell[q_point] = 0.0;
  for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
  {
  m_Nx_phasefield[q_point][k] = 0.0;
  m_grad_Nx_phasefield[q_point][k] = 0.0;
  m_Nx_disp[q_point][k] = 0.0;
  m_grad_Nx_disp[q_point][k] = 0.0;
  m_symm_grad_Nx_disp[q_point][k] = 0.0;
  }
  }
  }
  };
 
 
  template <int dim>
  struct PhaseFieldMonolithicSolve<dim>::ScratchData_ASM_RHS_BFGS
  {
  FEValues<dim> m_fe_values;
  FEFaceValues<dim> m_fe_face_values;
 
  std::vector<std::vector<double>> m_Nx_phasefield; // shape function values for phase-field
  std::vector<std::vector<Tensor<1, dim>>> m_grad_Nx_phasefield; // gradient of shape function values for phase field
 
  std::vector<std::vector<Tensor<1, dim>>> m_Nx_disp; // shape function values for displacement
  std::vector<std::vector<Tensor<2, dim>>> m_grad_Nx_disp; // gradient of shape function values for displacement
  std::vector<std::vector<SymmetricTensor<2, dim>>> m_symm_grad_Nx_disp; // symmetric gradient of shape function values for displacement
 
  const BlockVector<double>& m_solution_previous_step;
  std::vector<double> m_phasefield_previous_step_cell;
 
  ScratchData_ASM_RHS_BFGS(const FiniteElement<dim> & fe_cell,
  const QGauss<dim> & qf_cell,
  const UpdateFlags uf_cell,
  const QGauss<dim - 1> & qf_face,
  const UpdateFlags uf_face,
  const BlockVector<double>& solution_old)
  : m_fe_values(fe_cell, qf_cell, uf_cell)
  , m_fe_face_values(fe_cell, qf_face, uf_face)
  , m_Nx_phasefield(qf_cell.size(),
  std::vector<double>(fe_cell.n_dofs_per_cell()))
  , m_grad_Nx_phasefield(qf_cell.size(),
  std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
  , m_Nx_disp(qf_cell.size(),
  std::vector<Tensor<1, dim>>(fe_cell.n_dofs_per_cell()))
  , m_grad_Nx_disp(qf_cell.size(),
  std::vector<Tensor<2, dim>>(fe_cell.n_dofs_per_cell()))
  , m_symm_grad_Nx_disp(qf_cell.size(),
  std::vector<SymmetricTensor<2, dim>>(fe_cell.n_dofs_per_cell()))
  , m_solution_previous_step(solution_old)
  , m_phasefield_previous_step_cell(qf_cell.size())
  {}
 
  ScratchData_ASM_RHS_BFGS(const ScratchData_ASM_RHS_BFGS &rhs)
  : m_fe_values(rhs.m_fe_values.get_fe(),
  rhs.m_fe_values.get_quadrature(),
  rhs.m_fe_values.get_update_flags())
  , m_fe_face_values(rhs.m_fe_face_values.get_fe(),
  rhs.m_fe_face_values.get_quadrature(),
  rhs.m_fe_face_values.get_update_flags())
  , m_Nx_phasefield(rhs.m_Nx_phasefield)
  , m_grad_Nx_phasefield(rhs.m_grad_Nx_phasefield)
  , m_Nx_disp(rhs.m_Nx_disp)
  , m_grad_Nx_disp(rhs.m_grad_Nx_disp)
  , m_symm_grad_Nx_disp(rhs.m_symm_grad_Nx_disp)
  , m_solution_previous_step(rhs.m_solution_previous_step)
  , m_phasefield_previous_step_cell(rhs.m_phasefield_previous_step_cell)
  {}
 
  void reset()
  {
  const unsigned int n_q_points = m_Nx_phasefield.size();
  const unsigned int n_dofs_per_cell = m_Nx_phasefield[0].size();
  for (unsigned int q_point = 0; q_point < n_q_points; ++q_point)
  {
  Assert(m_Nx_phasefield[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_grad_Nx_phasefield[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  Assert(m_symm_grad_Nx_disp[q_point].size() == n_dofs_per_cell,
  ExcInternalError());
 
  m_phasefield_previous_step_cell[q_point] = 0.0;
  for (unsigned int k = 0; k < n_dofs_per_cell; ++k)
  {
  m_Nx_phasefield[q_point][k] = 0.0;
  m_grad_Nx_phasefield[q_point][k] = 0.0;
  m_Nx_disp[q_point][k] = 0.0;
  m_grad_Nx_disp[q_point][k] = 0.0;
  m_symm_grad_Nx_disp[q_point][k] = 0.0;
  }
  }
  }
  };
 
AssertThrow
#define AssertThrow(cond, exc)
Definition exception_macros.h:647
UpdateFlags
UpdateFlags
Definition fe_update_flags.h:62
update_values
@ update_values
Shape function values.
Definition fe_update_flags.h:73
update_gradients
@ update_gradients
Shape function gradients.
Definition fe_update_flags.h:79
size
std::size_t size
Definition mpi.cc:733
PlasticityLab::get_update_flags
*  *  *  *  *  *  TimeRateUpdateFlags TimeRateRequest< ValueType, dim, Number >  get_update_flags() const
Definition ALE_Finite_Strain_Plasticity.h:9929
PlasticityLab::const
*  *  *  *  std::vector< Number > ThermoPlasticMaterial< dim, ViscoplasticYieldLaw, Number >::get_state_parameters   const
Definition ALE_Finite_Strain_Plasticity.h:9082
WorkStream::run
void run(const std::vector< std::vector< Iterator > > &colored_iterators, Worker worker, Copier copier, const ScratchData &sample_scratch_data, const CopyData &sample_copy_data, const unsigned int queue_length=2 *MultithreadInfo::n_threads(), const unsigned int chunk_size=8)
Definition work_stream.h:1709
std
STL namespace.
types::material_id
unsigned int material_id
Definition types.h:182

constructor has no return type

  template <int dim>
  PhaseFieldMonolithicSolve<dim>::PhaseFieldMonolithicSolve(const std::string &input_file)
  : m_parameters(input_file)
  , m_triangulation(Triangulation<dim>::maximum_smoothing)
  , m_time(m_parameters.m_end_time)
  , m_logfile(m_parameters.m_logfile_name)
  , m_timer(m_logfile, TimerOutput::summary, TimerOutput::wall_times)
  , m_dof_handler(m_triangulation)
  , m_fe(FE_Q<dim>(m_parameters.m_poly_degree),
  dim, // displacement
  FE_Q<dim>(m_parameters.m_poly_degree),
  1) // phasefield
  , m_dofs_per_cell(m_fe.n_dofs_per_cell())
  , m_u_fe(m_first_u_component)
  , m_d_fe(m_d_component)
  , m_dofs_per_block(m_n_blocks)
  , m_qf_cell(m_parameters.m_quad_order)
  , m_qf_face(m_parameters.m_quad_order)
  , m_n_q_points(m_qf_cell.size())
  , m_vol_reference(0.0)
  {}
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid()
  {
  if (m_parameters.m_scenario == 1)
  make_grid_case_1();
  else if (m_parameters.m_scenario == 2)
  make_grid_case_2();
  else if (m_parameters.m_scenario == 3)
  make_grid_case_3();
  else if (m_parameters.m_scenario == 4)
  make_grid_case_4();
  else if (m_parameters.m_scenario == 5)
  make_grid_case_5();
  else if (m_parameters.m_scenario == 6)
  make_grid_case_6();
  else if (m_parameters.m_scenario == 7)
  make_grid_case_7();
  else if (m_parameters.m_scenario == 8)
  make_grid_case_8();
  else if (m_parameters.m_scenario == 9)
  make_grid_case_9();
  else if (m_parameters.m_scenario == 11)
  make_grid_case_11();
  else
  Assert(false, ExcMessage("The scenario has not been implemented!"));
 
  m_logfile << "\t\tTriangulation:"
  << "\n\t\t\tNumber of active cells: "
  << m_triangulation.n_active_cells()
  << "\n\t\t\tNumber of used vertices: "
  << m_triangulation.n_used_vertices()
  << std::endl;
 
  std::ofstream out("original_mesh.vtu");
  GridOut grid_out;
  grid_out.write_vtu(m_triangulation, out);
 
  m_vol_reference = GridTools::volume(m_triangulation);
  m_logfile << "\t\tGrid:\n\t\t\tReference volume: " << m_vol_reference << std::endl;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_1()
  {
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\tSquare tension (unstructured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("square_tension_unstructured.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] + 0.5 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 0.5 ) < 1.0e-9)
  face->set_boundary_id(1);
  else
  face->set_boundary_id(2);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[1]) < 0.01
  && cell->center()[0] > 0.495)
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[1] - 0.0) < 0.05
  && std::fabs(cell->center()[0] - 0.5) < 0.05)
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_2()
  {
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tSquare shear (unstructured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("square_shear_unstructured.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] + 0.5 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 0.5 ) < 1.0e-9)
  face->set_boundary_id(1);
  else if ( (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9)
  || (std::fabs(face->center()[0] - 1.0 ) < 1.0e-9))
  face->set_boundary_id(2);
  else
  face->set_boundary_id(3);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (cell->center()[0] > 0.45)
  && (cell->center()[1] < 0.05) )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 0.5) < 0.025
  && cell->center()[1] < 0.0 && cell->center()[1] > -0.025)
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_3()
  {
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\tSquare tension (structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("square_tension_structured.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 1.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else
  face->set_boundary_id(2);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (std::fabs(cell->center()[1] - 0.5) < 0.025)
  && (cell->center()[0] > 0.475) )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 0.5) < 0.025
  && std::fabs(cell->center()[1] - 0.5) < 0.025 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_4()
  {
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tSquare shear (structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("square_shear_structured.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 1.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else if ( (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9)
  || (std::fabs(face->center()[0] - 1.0 ) < 1.0e-9))
  face->set_boundary_id(2);
  else
  face->set_boundary_id(3);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (cell->center()[0] > 0.475)
  && (cell->center()[1] < 0.525) )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 0.5) < 0.025
  && cell->center()[1] < 0.5 && cell->center()[1] > 0.475 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_5()
  {
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tThree-point bending (structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("three_point_bending_structured.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 2.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else
  face->set_boundary_id(2);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 4.0) < 0.075
  && cell->center()[1] < 1.6)
  {
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
 
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 4.0) < 0.05
  && cell->center()[1] < 1.6)
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( std::fabs(cell->center()[0] - 4.0) < 0.075
  && std::fabs(cell->center()[1] - 0.4) < 0.075 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_6()
  {
  AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tSphere inclusion (3D structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  Triangulation<dim> tria_inner;
  GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.5);
 
  Triangulation<dim> tria_outer;
  GridGenerator::hyper_shell(
  tria_outer, Point<dim>(), 0.5, std::sqrt(dim), 2 * dim);
 
  Triangulation<dim> tmp_triangulation;
 
  GridGenerator::merge_triangulations(tria_inner, tria_outer, tmp_triangulation);
 
  tmp_triangulation.reset_all_manifolds();
  tmp_triangulation.set_all_manifold_ids(0);
 
  for (const auto &cell : tmp_triangulation.cell_iterators())
  {
  for (const auto &face : cell->face_iterators())
  {
  bool face_at_sphere_boundary = true;
  for (const auto v : face->vertex_indices())
  {
  if (std::abs(face->vertex(v).norm_square() - 0.25) > 1e-12)
  face_at_sphere_boundary = false;
  }
  if (face_at_sphere_boundary)
  face->set_all_manifold_ids(1);
  }
  if (cell->center().norm_square() < 0.25)
  cell->set_material_id(1);
  else
  cell->set_material_id(0);
  }
 
  tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 
  TransfiniteInterpolationManifold<dim> transfinite_manifold;
  transfinite_manifold.initialize(tmp_triangulation);
  tmp_triangulation.set_manifold(0, transfinite_manifold);
 
  tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  std::set<typename Triangulation< dim >::active_cell_iterator >
  cells_to_remove;
 
  for (const auto &cell : tmp_triangulation.active_cell_iterators())
  {
  if ( cell->center()[0] < 0.0
  || cell->center()[1] < 0.0
  || cell->center()[2] < 0.0)
  {
  cells_to_remove.insert(cell);
  }
  }
 
  GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
  cells_to_remove,
  m_triangulation);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(2);
  else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
  face->set_boundary_id(3);
  else
  face->set_boundary_id(4);
  }
  }
 
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( cell->center()[2] > 0.525
  && cell->center()[2] < 0.575
  && cell->center()[0] < 0.05
  && cell->center()[1] < 0.05 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::cbrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_7()
  {
  AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tSphere inclusion (3D structured version 2)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  Triangulation<dim> tria_inner;
  GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.49);
 
  Triangulation<dim> tria_outer;
  GridGenerator::hyper_shell(
  tria_outer, Point<dim>(), 0.49, std::sqrt(dim)*0.5, 2 * dim);
 
  Triangulation<dim> cube1;
  GridGenerator::hyper_rectangle(cube1, Point<dim>(0, 0, 0.5), Point<dim>(1, 1, 1.5));
  Triangulation<dim> cube2;
  GridGenerator::hyper_rectangle(cube2, Point<dim>(0, 0.5, -0.5), Point<dim>(1, 1.5, 0.5));
  Triangulation<dim> cube3;
  GridGenerator::hyper_rectangle(cube3, Point<dim>(0.5, -0.5, -0.5), Point<dim>(1.5, 0.5, 0.5));
 
  Triangulation<dim> tmp_triangulation;
  GridGenerator::merge_triangulations({&tria_inner, &tria_outer,
  &cube1, &cube2, &cube3}, tmp_triangulation);
 
  tmp_triangulation.reset_all_manifolds();
  tmp_triangulation.set_all_manifold_ids(0);
 
  for (const auto &cell : tmp_triangulation.cell_iterators())
  {
  for (const auto &face : cell->face_iterators())
  {
  bool face_at_sphere_boundary = true;
  for (const auto v : face->vertex_indices())
  {
  if (std::abs(face->vertex(v).norm_square() - 0.49 * 0.49) > 1e-12)
  face_at_sphere_boundary = false;
  }
  if (face_at_sphere_boundary)
  face->set_all_manifold_ids(1);
  }
  if (cell->center().norm_square() < 0.1)
  cell->set_material_id(1);
  else
  cell->set_material_id(0);
  }
 
  tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 
  TransfiniteInterpolationManifold<dim> transfinite_manifold;
  transfinite_manifold.initialize(tmp_triangulation);
  tmp_triangulation.set_manifold(0, transfinite_manifold);
 
  tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  std::set<typename Triangulation< dim >::active_cell_iterator >
  cells_to_remove;
 
  for (const auto &cell : tmp_triangulation.active_cell_iterators())
  {
  if ( cell->center()[0] < 0.0
  || cell->center()[1] < 0.0
  || cell->center()[2] < 0.0
  || cell->center()[0] > 1.0
  || cell->center()[1] > 1.0
  || cell->center()[2] > 1.0)
  {
  cells_to_remove.insert(cell);
  }
  }
 
  GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
  cells_to_remove,
  m_triangulation);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(2);
  else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
  face->set_boundary_id(3);
  else
  face->set_boundary_id(4);
  }
  }
 
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( cell->center()[2] > 0.505
  && cell->center()[2] < 0.575
  && cell->center()[0] < 0.05
  && cell->center()[1] < 0.05 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::cbrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_8()
  {
  AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tSphere inclusion (3D structured version 2 with barriers)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  Triangulation<dim> tria_inner;
  GridGenerator::hyper_ball(tria_inner, Point<dim>(), 0.49);
 
  Triangulation<dim> tria_outer;
  GridGenerator::hyper_shell(
  tria_outer, Point<dim>(), 0.49, std::sqrt(dim)*0.5, 2 * dim);
 
  Triangulation<dim> cube1;
  GridGenerator::hyper_rectangle(cube1, Point<dim>(0, 0, 0.5), Point<dim>(1, 1, 1.5));
  Triangulation<dim> cube2;
  GridGenerator::hyper_rectangle(cube2, Point<dim>(0, 0.5, -0.5), Point<dim>(1, 1.5, 0.5));
  Triangulation<dim> cube3;
  GridGenerator::hyper_rectangle(cube3, Point<dim>(0.5, -0.5, -0.5), Point<dim>(1.5, 0.5, 0.5));
 
  Triangulation<dim> tmp_triangulation;
  GridGenerator::merge_triangulations({&tria_inner, &tria_outer,
  &cube1, &cube2, &cube3}, tmp_triangulation);
 
  tmp_triangulation.reset_all_manifolds();
  tmp_triangulation.set_all_manifold_ids(0);
 
  for (const auto &cell : tmp_triangulation.cell_iterators())
  {
  for (const auto &face : cell->face_iterators())
  {
  bool face_at_sphere_boundary = true;
  for (const auto v : face->vertex_indices())
  {
  if (std::abs(face->vertex(v).norm_square() - 0.49 * 0.49) > 1e-12)
  face_at_sphere_boundary = false;
  }
  if (face_at_sphere_boundary)
  face->set_all_manifold_ids(1);
  }
  if (cell->center().norm_square() < 0.1)
  cell->set_material_id(1);
  else
  cell->set_material_id(0);
  }
 
  tmp_triangulation.set_manifold(1, SphericalManifold<dim>());
 
  TransfiniteInterpolationManifold<dim> transfinite_manifold;
  transfinite_manifold.initialize(tmp_triangulation);
  tmp_triangulation.set_manifold(0, transfinite_manifold);
 
  tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 
for
*  *  for(const auto &cell :triangulation.active_cell_iterators())
Definition ALE_Finite_Strain_Plasticity.h:8203
FE_Q
Definition fe_q.h:552
GridIn::attach_triangulation
void attach_triangulation(Triangulation< dim, spacedim > &tria)
Definition grid_in.cc:155
GridOut::write_vtu
void write_vtu(const Triangulation< dim, spacedim > &tria, std::ostream &out) const
Definition grid_out.cc:3557
TransfiniteInterpolationManifold::initialize
void initialize(const Triangulation< dim, spacedim > &triangulation)
Definition manifold_lib.cc:1656
vertex_indices
unsigned int vertex_indices[2]
Definition grid_tools_topology.cc:941
Triangulation::reset_all_manifolds
void reset_all_manifolds()
Differentiation::SD::fabs
Expression fabs(const Expression &x)
Definition symengine_math.cc:272
GridGenerator::hyper_shell
void hyper_shell(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center, const double inner_radius, const double outer_radius, const unsigned int n_cells=0, bool colorize=false)
GridGenerator::hyper_rectangle
void hyper_rectangle(Triangulation< dim, spacedim > &tria, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)
GridGenerator::hyper_ball
void hyper_ball(Triangulation< dim, spacedim > &tria, const Point< spacedim > &center={}, const double radius=1., const bool attach_spherical_manifold_on_boundary_cells=false)
GridGenerator::create_triangulation_with_removed_cells
void create_triangulation_with_removed_cells(const Triangulation< dim, spacedim > &input_triangulation, const std::set< typename Triangulation< dim, spacedim >::active_cell_iterator > &cells_to_remove, Triangulation< dim, spacedim > &result)
GridGenerator::merge_triangulations
void merge_triangulations(const Triangulation< dim, spacedim > &triangulation_1, const Triangulation< dim, spacedim > &triangulation_2, Triangulation< dim, spacedim > &result, const double duplicated_vertex_tolerance=1.0e-12, const bool copy_manifold_ids=false, const bool copy_boundary_ids=false)
GridTools::volume
double volume(const Triangulation< dim, spacedim > &tria)
Definition grid_tools_geometry.cc:98
Physics::Elasticity::Kinematics::e
SymmetricTensor< 2, dim, Number > e(const Tensor< 2, dim, Number > &F)
PlasticityLab::if
*  *  if(update_pressure &update_flags) *  compute_pressure(constitutive_request
std::sqrt
::VectorizedArray< Number, width > sqrt(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6505
std::abs
::VectorizedArray< Number, width > abs(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6563

some extra barriers

  for (const auto &cell : tmp_triangulation.cell_iterators())
  {
  if ( std::fabs(cell->center()[1] - 0.75) < 0.05
  && std::fabs(cell->center()[2] - 0.5625) < 0.05
  && std::fabs(cell->center()[0] - 0.0) < 0.2)
  cell->set_material_id(1);
 
  if ( std::fabs(cell->center()[1] - 0.0) < 0.2
  && std::fabs(cell->center()[2] - 0.5) < 0.1
  && std::fabs(cell->center()[0] - 0.75) < 0.05)
  cell->set_material_id(1);
  }
 
  std::set<typename Triangulation< dim >::active_cell_iterator >
  cells_to_remove;
 
  for (const auto &cell : tmp_triangulation.active_cell_iterators())
  {
  if ( cell->center()[0] < 0.0
  || cell->center()[1] < 0.0
  || cell->center()[2] < 0.0
  || cell->center()[0] > 1.0
  || cell->center()[1] > 1.0
  || cell->center()[2] > 1.0)
  {
  cells_to_remove.insert(cell);
  }
  }
 
  GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
  cells_to_remove,
  m_triangulation);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[0] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(1);
  else if (std::fabs(face->center()[2] - 0.0 ) < 1.0e-9)
  face->set_boundary_id(2);
  else if (std::fabs(face->center()[2] - 1.0 ) < 1.0e-9)
  face->set_boundary_id(3);
  else
  face->set_boundary_id(4);
  }
  }
 
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( cell->center()[2] > 0.505
  && cell->center()[2] < 0.575
  && cell->center()[0] < 0.05
  && cell->center()[1] < 0.05 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::cbrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_9()
  {
  AssertThrow(dim==2, ExcMessage("The dimension has to be 2D!"));
 
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tL-shape bending (2D structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  GridIn<dim> gridin;
  gridin.attach_triangulation(m_triangulation);
  std::ifstream f("L-Shape.msh");
  gridin.read_msh(f);
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[1] - 0.0 ) < 1.0e-9 )
  face->set_boundary_id(0);
  else
  face->set_boundary_id(1);
  }
  }
 
  m_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  unsigned int material_id;
  double length_scale;
  for (unsigned int i = 0; i < m_parameters.m_local_prerefine_times; i++)
  {
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (cell->center()[1] > 242.0)
  && (cell->center()[1] < 312.5)
  && (cell->center()[0] < 258.0) )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  cell->set_refine_flag();
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (cell->center()[0] - 250) < 0.0
  && (cell->center()[0] - 240) > 0.0
  && std::fabs(cell->center()[1] - 250) < 10.0 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::sqrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_grid_case_11()
  {
  AssertThrow(dim==3, ExcMessage("The dimension has to be 3D!"));
 
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
  m_logfile << "\t\t\t\tBrokenshire torsion (3D structured)" << std::endl;
  for (unsigned int i = 0; i < 80; ++i)
  m_logfile << "*";
  m_logfile << std::endl;
 
  Triangulation<2> triangulation_2d;
 
  double const length = 200.0;
  double const width = 50.0;
  double const height = 50.0;
  double const delta_L = 25.0;
  double const tan_theta = delta_L / (0.5*width);
 
  std::vector<unsigned int> repetitions(2, 1);
  repetitions[0] = 20;
  repetitions[1] = 5;
 
  Point<2> point1(0.0, 0.0);
  Point<2> point2(length, width);
 
  GridGenerator::subdivided_hyper_rectangle(triangulation_2d,
  repetitions,
  point1,
  point2 );
 
  typename Triangulation<2>::vertex_iterator vertex_ptr;
  vertex_ptr = triangulation_2d.begin_active_vertex();
  while (vertex_ptr != triangulation_2d.end_vertex())
  {
  Point<2> & vertex_point = vertex_ptr->vertex();
 
  const double delta_x = (vertex_point(1) - 0.5*width) * tan_theta;
 
  if (std::fabs(vertex_point(0) - 0.5*length) < 1.0e-6)
  {
  vertex_point(0) += delta_x;
  }
  else if (std::fabs(vertex_point(0) + length/repetitions[0] - 0.5*length) < 1.0e-6)
  {
  vertex_point(0) += (delta_x + length/repetitions[0]*0.5);
  }
  else if (std::fabs(vertex_point(0) - length/repetitions[0] - 0.5*length) < 1.0e-6)
  {
  vertex_point(0) += (delta_x - length/repetitions[0]*0.5);
  }
  else if (vertex_point(0) < 0.5*length - length/repetitions[0] - 1.0e-6)
  {
  vertex_point(0) += (delta_x + length/repetitions[0]*0.5) * vertex_point(0)/(0.5*length - length/repetitions[0]);
  }
  else if (vertex_point(0) > 0.5*length + length/repetitions[0] + 1.0e-6)
  {
  vertex_point(0) += (delta_x - length/repetitions[0]*0.5) * (length - vertex_point(0))/(0.5*length - length/repetitions[0]);
  }
 
  ++vertex_ptr;
  }
 
  Triangulation<dim> tmp_triangulation;
  const unsigned int n_layer = repetitions[1] + 1;
  GridGenerator::extrude_triangulation(triangulation_2d, n_layer, height, tmp_triangulation);
 
  tmp_triangulation.refine_global(m_parameters.m_global_refine_times);
 
  std::set<typename Triangulation< dim >::active_cell_iterator >
  cells_to_remove;
 
  for (const auto &cell : tmp_triangulation.active_cell_iterators())
  {
  if ( (std::fabs(cell->center()[0] - (cell->center()[1] - 0.5*width)*tan_theta - 0.5*length) < 2.5)
  && cell->center()[2] > 0.5* height )
  {
  cells_to_remove.insert(cell);
  }
  }
 
  GridGenerator::create_triangulation_with_removed_cells(tmp_triangulation,
  cells_to_remove,
  m_triangulation);
 
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
  unsigned int material_id;
  double length_scale;
  bool initiation_point_refine_unfinished = true;
  while (initiation_point_refine_unfinished)
  {
  initiation_point_refine_unfinished = false;
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  if ( (std::fabs(cell->center()[0] - (cell->center()[1] - 0.5*width)*tan_theta - 0.5*length) < 5.0)
  && cell->center()[2] <= 0.5*height
  && cell->center()[2] > 0.5*height - 5.0 )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if ( std::cbrt(cell->measure())
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  cell->set_refine_flag();
  initiation_point_refine_unfinished = true;
  }
  }
  }
  m_triangulation.execute_coarsening_and_refinement();
  }
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
 
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() == true)
  {
  if (std::fabs(face->center()[0] - length) < 1.0e-6 )
  face->set_boundary_id(0);
  else if (std::fabs(face->center()[0] - 0.0) < 1.0e-6 )
  face->set_boundary_id(1);
  else
  face->set_boundary_id(2);
  }
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::setup_system()
  {
  m_timer.enter_subsection("Setup system");
 
  std::vector<unsigned int> block_component(m_n_components,
  m_u_dof); // displacement
  block_component[m_d_component] = m_d_dof; // phasefield
 
  m_dof_handler.distribute_dofs(m_fe);
  DoFRenumbering::Cuthill_McKee(m_dof_handler);
  DoFRenumbering::component_wise(m_dof_handler, block_component);
 
  m_constraints.clear();
  DoFTools::make_hanging_node_constraints(m_dof_handler, m_constraints);
  m_constraints.close();
 
  m_dofs_per_block =
  DoFTools::count_dofs_per_fe_block(m_dof_handler, block_component);
 
  m_logfile << "\t\tTriangulation:"
  << "\n\t\t\t Number of active cells: "
  << m_triangulation.n_active_cells()
  << "\n\t\t\t Number of used vertices: "
  << m_triangulation.n_used_vertices()
  << "\n\t\t\t Number of active edges: "
  << m_triangulation.n_active_lines()
  << "\n\t\t\t Number of active faces: "
  << m_triangulation.n_active_faces()
  << "\n\t\t\t Number of degrees of freedom (total): "
  << m_dof_handler.n_dofs()
  << "\n\t\t\t Number of degrees of freedom (disp): "
  << m_dofs_per_block[m_u_dof]
  << "\n\t\t\t Number of degrees of freedom (phasefield): "
  << m_dofs_per_block[m_d_dof]
  << std::endl;
 
  m_tangent_matrix.clear();
  {
  BlockDynamicSparsityPattern dsp(m_dofs_per_block, m_dofs_per_block);
 
  Table<2, DoFTools::Coupling> coupling(m_n_components, m_n_components);
  for (unsigned int ii = 0; ii < m_n_components; ++ii)
  for (unsigned int jj = 0; jj < m_n_components; ++jj)
  coupling[ii][jj] = DoFTools::always;
 
  DoFTools::make_sparsity_pattern(
  m_dof_handler, coupling, dsp, m_constraints, false);
  m_sparsity_pattern.copy_from(dsp);
  }
 
  m_tangent_matrix.reinit(m_sparsity_pattern);
 
  m_system_rhs.reinit(m_dofs_per_block);
  m_solution.reinit(m_dofs_per_block);
 
  setup_qph();
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::make_constraints(const unsigned int it_nr)
  {
  const bool apply_dirichlet_bc = (it_nr == 0);
 
  if (it_nr > 1)
  {
  if (m_parameters.m_output_iteration_history)
  m_logfile << " --- " << std::flush;
  return;
  }
 
  if (m_parameters.m_output_iteration_history)
  m_logfile << " CST " << std::flush;
 
  if (apply_dirichlet_bc)
  {
  m_constraints.clear();
  DoFTools::make_hanging_node_constraints(m_dof_handler,
  m_constraints);
 
  const FEValuesExtractors::Scalar x_displacement(0);
  const FEValuesExtractors::Scalar y_displacement(1);
  const FEValuesExtractors::Scalar z_displacement(2);
 
  const FEValuesExtractors::Vector displacements(0);
 
  if ( m_parameters.m_scenario == 1
  || m_parameters.m_scenario == 3)
  {
height
*  const Number height
Definition ALE_Finite_Strain_Plasticity.h:8297
DoFTools::make_hanging_node_constraints
void make_hanging_node_constraints(const DoFHandler< dim, spacedim > &dof_handler, AffineConstraints< number > &constraints)
Definition dof_tools_constraints.cc:3094
DoFTools::make_sparsity_pattern
void make_sparsity_pattern(const DoFHandler< dim, spacedim > &dof_handler, SparsityPatternBase &sparsity_pattern, const AffineConstraints< number > &constraints={}, const bool keep_constrained_dofs=true, const types::subdomain_id subdomain_id=numbers::invalid_subdomain_id)
Definition dof_tools_sparsity.cc:55
DoFRenumbering::component_wise
void component_wise(DoFHandler< dim, spacedim > &dof_handler, const std::vector< unsigned int > &target_component=std::vector< unsigned int >())
Definition dof_renumbering.cc:670
DoFRenumbering::Cuthill_McKee
void Cuthill_McKee(DoFHandler< dim, spacedim > &dof_handler, const bool reversed_numbering=false, const bool use_constraints=false, const std::vector< types::global_dof_index > &starting_indices=std::vector< types::global_dof_index >())
Definition dof_renumbering.cc:370
DoFTools::count_dofs_per_fe_block
std::vector< types::global_dof_index > count_dofs_per_fe_block(const DoFHandler< dim, spacedim > &dof, const std::vector< unsigned int > &target_block=std::vector< unsigned int >())
Definition dof_tools.cc:2366
GridGenerator::extrude_triangulation
void extrude_triangulation(const Triangulation< 2, 2 > &input, const unsigned int n_slices, const double height, Triangulation< 3, 3 > &result, const bool copy_manifold_ids=false, const std::vector< types::manifold_id > &manifold_priorities={})
GridGenerator::subdivided_hyper_rectangle
void subdivided_hyper_rectangle(Triangulation< dim, spacedim > &tria, const std::vector< unsigned int > &repetitions, const Point< dim > &p1, const Point< dim > &p2, const bool colorize=false)

Dirichlet B,C. bottom surface

  const int boundary_id_bottom_surface = 0;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_bottom_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(y_displacement));
 
  typename Triangulation<dim>::active_vertex_iterator vertex_itr;
  vertex_itr = m_triangulation.begin_active_vertex();
  std::vector<types::global_dof_index> node_xy(m_fe.dofs_per_vertex);
 
  for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
  {
  if ( (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
  && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
  {
  node_xy = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
  }
  }
  m_constraints.add_line(node_xy[0]);
  m_constraints.set_inhomogeneity(node_xy[0], 0.0);
 
  m_constraints.add_line(node_xy[1]);
  m_constraints.set_inhomogeneity(node_xy[1], 0.0);
 
  const int boundary_id_top_surface = 1;
  /*
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_top_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(x_displacement));
  */
  const double time_inc = m_time.get_delta_t();
  double disp_magnitude = m_time.get_magnitude();
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_top_surface,
  Functions::ConstantFunction<dim>(
  disp_magnitude*time_inc, m_n_components),
  m_constraints,
  m_fe.component_mask(y_displacement));
  }
  else if ( m_parameters.m_scenario == 2
  || m_parameters.m_scenario == 4)
  {
VectorTools::interpolate_boundary_values
void interpolate_boundary_values(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const std::map< types::boundary_id, const Function< spacedim, number > * > &function_map, std::map< types::global_dof_index, number > &boundary_values, const ComponentMask &component_mask={})

Dirichlet B,C. bottom surface

  const int boundary_id_bottom_surface = 0;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_bottom_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(displacements));
 
  const int boundary_id_top_surface = 1;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_top_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(y_displacement));
 
  const double time_inc = m_time.get_delta_t();
  double disp_magnitude = m_time.get_magnitude();
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_top_surface,
  Functions::ConstantFunction<dim>(
  disp_magnitude*time_inc, m_n_components),
  m_constraints,
  m_fe.component_mask(x_displacement));
 
  const int boundary_id_side_surfaces = 2;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_side_surfaces,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(y_displacement));
  }
  else if (m_parameters.m_scenario == 5)
  {
  typename Triangulation<dim>::active_vertex_iterator vertex_itr;
  vertex_itr = m_triangulation.begin_active_vertex();
  std::vector<types::global_dof_index> node_bottomleft(m_fe.dofs_per_vertex);
  std::vector<types::global_dof_index> node_bottomright(m_fe.dofs_per_vertex);
  std::vector<types::global_dof_index> node_topcenter(m_fe.dofs_per_vertex);
 
  for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
  {
  if ( (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
  && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
  {
  node_bottomleft = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
  }
  if ( (std::fabs(vertex_itr->vertex()[0] - 8.0) < 1.0e-9)
  && (std::fabs(vertex_itr->vertex()[1] - 0.0) < 1.0e-9) )
  {
  node_bottomright = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
  }
  if ( (std::fabs(vertex_itr->vertex()[0] - 4.0) < 1.0e-9)
  && (std::fabs(vertex_itr->vertex()[1] - 2.0) < 1.0e-9) )
  {
  node_topcenter = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
  }
  }

bottom-left node fixed in both x- and y-directions

  m_constraints.add_line(node_bottomleft[0]);
  m_constraints.set_inhomogeneity(node_bottomleft[0], 0.0);
 
  m_constraints.add_line(node_bottomleft[1]);
  m_constraints.set_inhomogeneity(node_bottomleft[1], 0.0);
 

bottom-right node only fixed in y-direction

  m_constraints.add_line(node_bottomright[1]);
  m_constraints.set_inhomogeneity(node_bottomright[1], 0.0);
 

top-center node applied with y-displacement

  const double time_inc = m_time.get_delta_t();
  double disp_magnitude = m_time.get_magnitude();
 
  m_constraints.add_line(node_topcenter[1]);
  m_constraints.set_inhomogeneity(node_topcenter[1], disp_magnitude*time_inc);
  }
  else if ( m_parameters.m_scenario == 6
  || m_parameters.m_scenario == 7
  || m_parameters.m_scenario == 8)
  {
  const int x0_surface = 0;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  x0_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(x_displacement));
  const int y0_surface = 1;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  y0_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(y_displacement));
  const int z0_surface = 2;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  z0_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(z_displacement));
 
  const int z1_surface = 3;
  const double time_inc = m_time.get_delta_t();
  double disp_magnitude = m_time.get_magnitude();
  VectorTools::interpolate_boundary_values(m_dof_handler,
  z1_surface,
  Functions::ConstantFunction<dim>(
  disp_magnitude*time_inc, m_n_components),
  m_constraints,
  m_fe.component_mask(z_displacement));
  }
  else if (m_parameters.m_scenario == 9)
  {

Dirichlet B,C. bottom surface

  const int boundary_id_bottom_surface = 0;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_bottom_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(displacements));
 
  typename Triangulation<dim>::active_vertex_iterator vertex_itr;
  vertex_itr = m_triangulation.begin_active_vertex();
  std::vector<types::global_dof_index> node_disp_control(m_fe.dofs_per_vertex);
 
  for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
  {
  if ( (std::fabs(vertex_itr->vertex()[0] - 470.0) < 1.0e-9)
  && (std::fabs(vertex_itr->vertex()[1] - 250.0) < 1.0e-9) )
  {
  node_disp_control = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);

node applied with y-displacement

  const double time_inc = m_time.get_delta_t();
  double disp_magnitude = m_time.get_magnitude();
 
  m_constraints.add_line(node_disp_control[1]);
  m_constraints.set_inhomogeneity(node_disp_control[1], disp_magnitude*time_inc);
  }
  }
  }
  else if (m_parameters.m_scenario == 11)
  {

Dirichlet B,C. right surface

  const int boundary_id_right_surface = 0;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_right_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(displacements));
 

Dirichlet B,C. left surface

  const int boundary_id_left_surface = 1;
  VectorTools::interpolate_boundary_values(m_dof_handler,
  boundary_id_left_surface,
  Functions::ZeroFunction<dim>(m_n_components),
  m_constraints,
  m_fe.component_mask(x_displacement));
 
  typename Triangulation<dim>::active_vertex_iterator vertex_itr;
  vertex_itr = m_triangulation.begin_active_vertex();
  std::vector<types::global_dof_index> node_rotate(m_fe.dofs_per_vertex);
  double node_dist = 0.0;
  double disp_mag = 0.0;
  double angle_theta = 0.0;
  double disp_y = 0;
  double disp_z = 0;
 
  for (; vertex_itr != m_triangulation.end_vertex(); ++vertex_itr)
  {
  if (std::fabs(vertex_itr->vertex()[0] - 0.0) < 1.0e-9)
  {
  node_rotate = usr_utilities::get_vertex_dofs(vertex_itr, m_dof_handler);
  node_dist = std::sqrt( vertex_itr->vertex()[1] * vertex_itr->vertex()[1]
  + vertex_itr->vertex()[2] * vertex_itr->vertex()[2]);
 
  angle_theta = m_time.get_delta_t() * m_time.get_magnitude();
  disp_mag = node_dist * std::tan(angle_theta);
 
  if (node_dist > 0)
  {
  disp_y = vertex_itr->vertex()[2]/node_dist * disp_mag;
  disp_z = -vertex_itr->vertex()[1]/node_dist * disp_mag;
  }
  else
  {
  disp_y = 0.0;
  disp_z = 0.0;
  }
 
  m_constraints.add_line(node_rotate[1]);
  m_constraints.set_inhomogeneity(node_rotate[1], disp_y);
 
  m_constraints.add_line(node_rotate[2]);
  m_constraints.set_inhomogeneity(node_rotate[2], disp_z);
  }
  }
  }
  else
  Assert(false, ExcMessage("The scenario has not been implemented!"));
  }
  else // inhomogeneous constraints
  {
  if (m_constraints.has_inhomogeneities())
  {
  AffineConstraints<double> homogeneous_constraints(m_constraints);
  for (unsigned int dof = 0; dof != m_dof_handler.n_dofs(); ++dof)
  if (homogeneous_constraints.is_inhomogeneously_constrained(dof))
  homogeneous_constraints.set_inhomogeneity(dof, 0.0);
  m_constraints.clear();
  m_constraints.copy_from(homogeneous_constraints);
  }
  }
  m_constraints.close();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_newton(const BlockVector<double> & solution_old)
  {
  m_timer.enter_subsection("Assemble system");
 
  if (m_parameters.m_output_iteration_history)
  m_logfile << " ASM_SYS " << std::flush;
 
  m_tangent_matrix = 0.0;
  m_system_rhs = 0.0;
 
  const UpdateFlags uf_cell(update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const UpdateFlags uf_face(update_values | update_normal_vectors |
  update_JxW_values);
 
  PerTaskData_ASM per_task_data(m_fe.n_dofs_per_cell());
  ScratchData_ASM scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 
  auto worker =
  [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data)
  {
  this->assemble_system_newton_one_cell(cell, scratch, data);
  };
 
  auto copier = [this](const PerTaskData_ASM &data)
  {
  this->m_constraints.distribute_local_to_global(data.m_cell_matrix,
  data.m_cell_rhs,
  data.m_local_dof_indices,
  m_tangent_matrix,
  m_system_rhs);
  };
 
  WorkStream::run(
  m_dof_handler.active_cell_iterators(),
  worker,
  copier,
  scratch_data,
  per_task_data);
 
  m_timer.leave_subsection();
  }
 
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_B0(const BlockVector<double> & solution_old)
  {
  m_timer.enter_subsection("Assemble B0");
 
  m_tangent_matrix = 0.0;
 
  const UpdateFlags uf_cell(update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const UpdateFlags uf_face(update_values | update_normal_vectors |
  update_JxW_values);
 
  PerTaskData_ASM per_task_data(m_fe.n_dofs_per_cell());
  ScratchData_ASM scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 
  auto worker =
  [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data)
  {
  this->assemble_system_B0_one_cell(cell, scratch, data);
  };
 
  auto copier = [this](const PerTaskData_ASM &data)
  {
  this->m_constraints.distribute_local_to_global(data.m_cell_matrix,
  data.m_local_dof_indices,
  m_tangent_matrix);
  };
 
  WorkStream::run(
  m_dof_handler.active_cell_iterators(),
  worker,
  copier,
  scratch_data,
  per_task_data);
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS_parallel(const BlockVector<double> & solution_old,
  BlockVector<double> & system_rhs)
  {
  m_timer.enter_subsection("Assemble RHS");
 
update_normal_vectors
@ update_normal_vectors
Normal vectors.
Definition fe_update_flags.h:139
update_JxW_values
@ update_JxW_values
Transformed quadrature weights.
Definition fe_update_flags.h:132
update_quadrature_points
@ update_quadrature_points
Transformed quadrature points.
Definition fe_update_flags.h:125
std::tan
::VectorizedArray< Number, width > tan(const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6265

m_logfile << " A_RHS " << std::flush;

  system_rhs = 0.0;
 
  const UpdateFlags uf_cell(update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const UpdateFlags uf_face(update_values | update_normal_vectors |
  update_JxW_values);
 
  PerTaskData_ASM_RHS_BFGS per_task_data(m_fe.n_dofs_per_cell());
  ScratchData_ASM_RHS_BFGS scratch_data(m_fe, m_qf_cell, uf_cell, m_qf_face, uf_face, solution_old);
 
  auto worker =
  [this](const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM_RHS_BFGS & scratch,
  PerTaskData_ASM_RHS_BFGS & data)
  {
  this->assemble_system_rhs_BFGS_one_cell(cell, scratch, data);
  };
 
  auto copier = [this, &system_rhs](const PerTaskData_ASM_RHS_BFGS &data)
  {
  this->m_constraints.distribute_local_to_global(data.m_cell_rhs,
  data.m_local_dof_indices,
  system_rhs);
  };
 
  WorkStream::run(
  m_dof_handler.active_cell_iterators(),
  worker,
  copier,
  scratch_data,
  per_task_data);
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM_RHS_BFGS & scratch,
  PerTaskData_ASM_RHS_BFGS & data) const
  {
  data.reset();
  scratch.reset();
  scratch.m_fe_values.reinit(cell);
  cell->get_dof_indices(data.m_local_dof_indices);
 
  scratch.m_fe_values[m_d_fe].get_function_values(
  scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 
  const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  const double time_ramp = (m_time.current() / m_time.end());
  std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 
  right_hand_side(scratch.m_fe_values.get_quadrature_points(),
  rhs_values,
  m_parameters.m_x_component*1.0,
  m_parameters.m_y_component*1.0,
  m_parameters.m_z_component*1.0);
 
  const double delta_time = m_time.get_delta_t();
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  for (const unsigned int k : scratch.m_fe_values.dof_indices())
  {
  const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 
  if (k_group == m_u_dof)
  {
  scratch.m_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].value(k, q_point);
  scratch.m_grad_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].gradient(k, q_point);
  scratch.m_symm_grad_Nx_disp[q_point][k] =
  symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
  }
  else if (k_group == m_d_dof)
  {
  scratch.m_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].value(k, q_point);
  scratch.m_grad_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].gradient(k, q_point);
  }
  else
  Assert(k_group <= m_d_dof, ExcInternalError());
  }
  }
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  const double length_scale = lqph[q_point]->get_length_scale();
  const double gc = lqph[q_point]->get_critical_energy_release_rate();
  const double eta = lqph[q_point]->get_viscosity();
  const double history_strain_energy = lqph[q_point]->get_history_max_positive_strain_energy();
  const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 
  double history_value = history_strain_energy;
  if (current_positive_strain_energy > history_strain_energy)
  history_value = current_positive_strain_energy;
 
  const double phasefield_value = lqph[q_point]->get_phase_field_value();
  const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 
  const std::vector<double> & N_phasefield = scratch.m_Nx_phasefield[q_point];
  const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
  const double old_phasefield = scratch.m_phasefield_previous_step_cell[q_point];
 
  const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 
  const std::vector<Tensor<1,dim>> & N_disp = scratch.m_Nx_disp[q_point];
  const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
  scratch.m_symm_grad_Nx_disp[q_point];
  const double JxW = scratch.m_fe_values.JxW(q_point);
 
  SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 
  for (const unsigned int i : scratch.m_fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  data.m_cell_rhs(i) += (symm_grad_N_disp[i] * cauchy_stress) * JxW;
 
first
Point< 2 > first
Definition grid_out.cc:4639

contributions from the body force to right-hand side

  data.m_cell_rhs(i) -= N_disp[i] * rhs_values[q_point] * JxW;
  }
  else if (i_group == m_d_dof)
  {
  data.m_cell_rhs(i) += ( gc * length_scale * grad_N_phasefield[i] * phasefield_grad
  + ( gc / length_scale * phasefield_value
  + eta / delta_time * (phasefield_value - old_phasefield)
  + degradation_function_derivative(phasefield_value) * history_value )
  * N_phasefield[i]
  ) * JxW;
  }
  else
  Assert(i_group <= m_d_dof, ExcInternalError());
  } // i
  } // q_point
 

if there is surface pressure, this surface pressure always applied to the reference configuration

  const unsigned int face_pressure_id = 100;
  const double p0 = 0.0;
 
  for (const auto &face : cell->face_iterators())
  if (face->at_boundary() && face->boundary_id() == face_pressure_id)
  {
  scratch.m_fe_face_values.reinit(cell, face);
 
  for (const unsigned int f_q_point : scratch.m_fe_face_values.quadrature_point_indices())
  {
  const Tensor<1, dim> &N = scratch.m_fe_face_values.normal_vector(f_q_point);
 
  const double pressure = p0 * time_ramp;
  const Tensor<1, dim> traction = pressure * N;
 
  for (const unsigned int i : scratch.m_fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  const unsigned int component_i = m_fe.system_to_component_index(i).first;
  const double Ni = scratch.m_fe_face_values.shape_value(i, f_q_point);
  const double JxW = scratch.m_fe_face_values.JxW(f_q_point);
  data.m_cell_rhs(i) -= (Ni * traction[component_i]) * JxW;
  }
  }
  }
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_newton_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data) const
  {
  data.reset();
  scratch.reset();
  scratch.m_fe_values.reinit(cell);
  cell->get_dof_indices(data.m_local_dof_indices);
 
  scratch.m_fe_values[m_d_fe].get_function_values(
  scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 
  const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  const double time_ramp = (m_time.current() / m_time.end());
  std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
 
  right_hand_side(scratch.m_fe_values.get_quadrature_points(),
  rhs_values,
  m_parameters.m_x_component*1.0,
  m_parameters.m_y_component*1.0,
  m_parameters.m_z_component*1.0);
 
  const double delta_time = m_time.get_delta_t();
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  for (const unsigned int k : scratch.m_fe_values.dof_indices())
  {
  const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 
  if (k_group == m_u_dof)
  {
  scratch.m_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].value(k, q_point);
  scratch.m_grad_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].gradient(k, q_point);
  scratch.m_symm_grad_Nx_disp[q_point][k] =
  symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
  }
  else if (k_group == m_d_dof)
  {
  scratch.m_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].value(k, q_point);
  scratch.m_grad_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].gradient(k, q_point);
  }
  else
  Assert(k_group <= m_d_dof, ExcInternalError());
  }
  }
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  const double length_scale = lqph[q_point]->get_length_scale();
  const double gc = lqph[q_point]->get_critical_energy_release_rate();
  const double eta = lqph[q_point]->get_viscosity();
  const double history_strain_energy = lqph[q_point]->get_history_max_positive_strain_energy();
  const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 
  double history_value = history_strain_energy;
  if (current_positive_strain_energy > history_strain_energy)
  history_value = current_positive_strain_energy;
 
  const double phasefield_value = lqph[q_point]->get_phase_field_value();
  const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 
  const std::vector<double> & N_phasefield = scratch.m_Nx_phasefield[q_point];
  const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
  const double old_phasefield = scratch.m_phasefield_previous_step_cell[q_point];
 
  const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
  const SymmetricTensor<2, dim> & cauchy_stress_positive = lqph[q_point]->get_cauchy_stress_positive();
  const SymmetricTensor<4, dim> & mechanical_C = lqph[q_point]->get_mechanical_C();
 
  const std::vector<Tensor<1,dim>> & N_disp = scratch.m_Nx_disp[q_point];
  const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
  scratch.m_symm_grad_Nx_disp[q_point];
  const double JxW = scratch.m_fe_values.JxW(q_point);
 
  SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 
  for (const unsigned int i : scratch.m_fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  data.m_cell_rhs(i) -= (symm_grad_N_disp[i] * cauchy_stress) * JxW;
 

contributions from the body force to right-hand side

  data.m_cell_rhs(i) += N_disp[i] * rhs_values[q_point] * JxW;
  }
  else if (i_group == m_d_dof)
  {
  data.m_cell_rhs(i) -= ( gc * length_scale * grad_N_phasefield[i] * phasefield_grad
  + ( gc / length_scale * phasefield_value
  + eta / delta_time * (phasefield_value - old_phasefield)
  + degradation_function_derivative(phasefield_value) * history_value )
  * N_phasefield[i]
  ) * JxW;
  }
  else
  Assert(i_group <= m_d_dof, ExcInternalError());
 
  if (i_group == m_u_dof)
  {
  symm_grad_Nx_i_x_C = symm_grad_N_disp[i] * mechanical_C;
  }
 
  for (const unsigned int j : scratch.m_fe_values.dof_indices())
  {
  const unsigned int j_group = m_fe.system_to_base_index(j).first.first;
 
  if ((i_group == j_group) && (i_group == m_u_dof))
  {
  data.m_cell_matrix(i, j) += symm_grad_Nx_i_x_C * symm_grad_N_disp[j] * JxW;
  }
  else if ((i_group == j_group) && (i_group == m_d_dof))
  {
  data.m_cell_matrix(i, j) += ( ( gc/length_scale + eta/delta_time +
  degradation_function_2nd_order_derivative(phasefield_value)
  * history_value )
  * N_phasefield[i] * N_phasefield[j]
  + gc * length_scale * grad_N_phasefield[i] * grad_N_phasefield[j]
  ) * JxW;
  }
  else if ((i_group == m_u_dof) && (j_group == m_d_dof))
  {
  data.m_cell_matrix(i, j) += symm_grad_N_disp[i] * cauchy_stress_positive
  * degradation_function_derivative(phasefield_value)
  * N_phasefield[j] * JxW;
  }
  else if ((i_group == m_d_dof) && (j_group == m_u_dof))
  {
  if (current_positive_strain_energy > history_strain_energy)
  data.m_cell_matrix(i, j) += N_phasefield[i]
  * degradation_function_derivative(phasefield_value)
  * cauchy_stress_positive
  * symm_grad_N_disp[j]
  * JxW;
  else
  data.m_cell_matrix(i, j) += 0.0;
  }
  else
  Assert((i_group <= m_d_dof) && (j_group <= m_d_dof),
  ExcInternalError());
  } // j
  } // i
  } // q_point
 

if there is surface pressure, this surface pressure always applied to the reference configuration

  const unsigned int face_pressure_id = 100;
  const double p0 = 0.0;
 
  for (const auto &face : cell->face_iterators())
  if (face->at_boundary() && face->boundary_id() == face_pressure_id)
  {
  scratch.m_fe_face_values.reinit(cell, face);
 
  for (const unsigned int f_q_point : scratch.m_fe_face_values.quadrature_point_indices())
  {
  const Tensor<1, dim> &N = scratch.m_fe_face_values.normal_vector(f_q_point);
 
  const double pressure = p0 * time_ramp;
  const Tensor<1, dim> traction = pressure * N;
 
  for (const unsigned int i : scratch.m_fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  const unsigned int component_i = m_fe.system_to_component_index(i).first;
  const double Ni = scratch.m_fe_face_values.shape_value(i, f_q_point);
  const double JxW = scratch.m_fe_face_values.JxW(f_q_point);
  data.m_cell_rhs(i) += (Ni * traction[component_i]) * JxW;
  }
  }
  }
  }
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_B0_one_cell(
  const typename DoFHandler<dim>::active_cell_iterator &cell,
  ScratchData_ASM & scratch,
  PerTaskData_ASM & data) const
  {
  data.reset();
  scratch.reset();
  scratch.m_fe_values.reinit(cell);
  cell->get_dof_indices(data.m_local_dof_indices);
 
  scratch.m_fe_values[m_d_fe].get_function_values(
  scratch.m_solution_previous_step, scratch.m_phasefield_previous_step_cell);
 
  const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  const double delta_time = m_time.get_delta_t();
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  for (const unsigned int k : scratch.m_fe_values.dof_indices())
  {
  const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 
  if (k_group == m_u_dof)
  {
  scratch.m_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].value(k, q_point);
  scratch.m_grad_Nx_disp[q_point][k] =
  scratch.m_fe_values[m_u_fe].gradient(k, q_point);
  scratch.m_symm_grad_Nx_disp[q_point][k] =
  symmetrize(scratch.m_grad_Nx_disp[q_point][k]);
  }
  else if (k_group == m_d_dof)
  {
  scratch.m_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].value(k, q_point);
  scratch.m_grad_Nx_phasefield[q_point][k] =
  scratch.m_fe_values[m_d_fe].gradient(k, q_point);
  }
  else
  Assert(k_group <= m_d_dof, ExcInternalError());
  }
  }
 
  for (const unsigned int q_point : scratch.m_fe_values.quadrature_point_indices())
  {
  const double length_scale = lqph[q_point]->get_length_scale();
  const double gc = lqph[q_point]->get_critical_energy_release_rate();
  const double eta = lqph[q_point]->get_viscosity();
  const double history_strain_energy = lqph[q_point]->get_history_max_positive_strain_energy();
  const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 
  double history_value = history_strain_energy;
  if (current_positive_strain_energy > history_strain_energy)
  history_value = current_positive_strain_energy;
 
  const double phasefield_value = lqph[q_point]->get_phase_field_value();
 
  const std::vector<double> & N_phasefield = scratch.m_Nx_phasefield[q_point];
  const std::vector<Tensor<1, dim>> & grad_N_phasefield = scratch.m_grad_Nx_phasefield[q_point];
 

const SymmetricTensor<2, dim> & cauchy_stress_positive = lqph[q_point]->get_cauchy_stress_positive();

  const SymmetricTensor<4, dim> & mechanical_C = lqph[q_point]->get_mechanical_C();
 
  const std::vector<SymmetricTensor<2, dim>> & symm_grad_N_disp =
  scratch.m_symm_grad_Nx_disp[q_point];
  const double JxW = scratch.m_fe_values.JxW(q_point);
 
  SymmetricTensor<2, dim> symm_grad_Nx_i_x_C;
 
  for (const unsigned int i : scratch.m_fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  symm_grad_Nx_i_x_C = symm_grad_N_disp[i] * mechanical_C;
  }
 
  for (const unsigned int j : scratch.m_fe_values.dof_indices())
  {
  const unsigned int j_group = m_fe.system_to_base_index(j).first.first;
 
  if ((i_group == j_group) && (i_group == m_u_dof))
  {
  data.m_cell_matrix(i, j) += symm_grad_Nx_i_x_C * symm_grad_N_disp[j] * JxW;
  }
  else if ((i_group == j_group) && (i_group == m_d_dof))
  {
  data.m_cell_matrix(i, j) += ( ( gc/length_scale + eta/delta_time
  + degradation_function_2nd_order_derivative(phasefield_value)
  * history_value )
  * N_phasefield[i] * N_phasefield[j]
  + gc * length_scale * grad_N_phasefield[i] * grad_N_phasefield[j]
  ) * JxW;
  }
  else
  Assert((i_group <= m_d_dof) && (j_group <= m_d_dof),
  ExcInternalError());
  } // j
  } // i
  } // q_point
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::assemble_system_rhs_BFGS(const BlockVector<double> & solution_old,
  BlockVector<double> & system_rhs)
  {
  m_timer.enter_subsection("Assemble RHS");
 

m_logfile << " A_RHS " << std::flush;

  system_rhs = 0.0;
 
  Vector<double> cell_rhs(m_dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices(m_dofs_per_cell);
 
  const double time_ramp = (m_time.current() / m_time.end());
  const double delta_time = m_time.get_delta_t();
 
  std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
  const UpdateFlags uf_cell(update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const UpdateFlags uf_face(update_values | update_normal_vectors |
  update_JxW_values);
 
  FEValues<dim> fe_values(m_fe, m_qf_cell, uf_cell);
  FEFaceValues<dim> fe_face_values(m_fe, m_qf_face, uf_face);
 

shape function values for displacement field

  std::vector<std::vector<Tensor<1, dim>>>
  Nx_disp(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
  std::vector<std::vector<Tensor<2, dim>>>
  grad_Nx_disp(m_qf_cell.size(), std::vector<Tensor<2, dim>>(m_dofs_per_cell));
  std::vector<std::vector<SymmetricTensor<2, dim>>>
  symm_grad_Nx_disp(m_qf_cell.size(), std::vector<SymmetricTensor<2, dim>>(m_dofs_per_cell));
 

shape function values for phase field

  std::vector<std::vector<double>>
  Nx_phasefield(m_qf_cell.size(), std::vector<double>(m_dofs_per_cell));
  std::vector<std::vector<Tensor<1, dim>>>
  grad_Nx_phasefield(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
 
  std::vector<double> phasefield_previous_step_cell(m_qf_cell.size());
 
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {
  const std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  cell_rhs = 0.0;
  fe_values.reinit(cell);
  right_hand_side(fe_values.get_quadrature_points(),
  rhs_values,
  m_parameters.m_x_component*time_ramp,
  m_parameters.m_y_component*time_ramp,
  m_parameters.m_z_component*time_ramp);
 
  fe_values[m_d_fe].get_function_values(
  solution_old, phasefield_previous_step_cell);
 
  for (const unsigned int q_point : fe_values.quadrature_point_indices())
  {
  for (const unsigned int k : fe_values.dof_indices())
  {
  const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 
  if (k_group == m_u_dof)
  {
  Nx_disp[q_point][k] = fe_values[m_u_fe].value(k, q_point);
  grad_Nx_disp[q_point][k] = fe_values[m_u_fe].gradient(k, q_point);
  symm_grad_Nx_disp[q_point][k] = symmetrize(grad_Nx_disp[q_point][k]);
  }
  else if (k_group == m_d_dof)
  {
  Nx_phasefield[q_point][k] = fe_values[m_d_fe].value(k, q_point);
  grad_Nx_phasefield[q_point][k] = fe_values[m_d_fe].gradient(k, q_point);
  }
  else
  Assert(k_group <= m_d_dof, ExcInternalError());
  }
  }
 
  for (const unsigned int q_point : fe_values.quadrature_point_indices())
  {
  const double length_scale = lqph[q_point]->get_length_scale();
  const double gc = lqph[q_point]->get_critical_energy_release_rate();
  const double eta = lqph[q_point]->get_viscosity();
  const double history_strain_energy = lqph[q_point]->get_history_max_positive_strain_energy();
  const double current_positive_strain_energy = lqph[q_point]->get_current_positive_strain_energy();
 
  double history_value = history_strain_energy;
  if (current_positive_strain_energy > history_strain_energy)
  history_value = current_positive_strain_energy;
 
  const double phasefield_value = lqph[q_point]->get_phase_field_value();
  const Tensor<1, dim> phasefield_grad = lqph[q_point]->get_phase_field_gradient();
 
  const std::vector<double> & N_phasefield = Nx_phasefield[q_point];
  const std::vector<Tensor<1, dim>> & grad_N_phasefield = grad_Nx_phasefield[q_point];
  const double old_phasefield = phasefield_previous_step_cell[q_point];
 
  const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 
  const std::vector<Tensor<1,dim>> & N = Nx_disp[q_point];
  const std::vector<SymmetricTensor<2, dim>> & symm_grad_N = symm_grad_Nx_disp[q_point];
  const double JxW = fe_values.JxW(q_point);
 
  for (const unsigned int i : fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  cell_rhs(i) += (symm_grad_N[i] * cauchy_stress) * JxW;
LAPACKSupport::N
constexpr char N
Definition lapack_support.h:155

contributions from the body force to right-hand side

  cell_rhs(i) -= N[i] * rhs_values[q_point] * JxW;
  }
  else if (i_group == m_d_dof)
  {
  cell_rhs(i) += ( gc * length_scale * grad_N_phasefield[i] * phasefield_grad
  + ( gc / length_scale * phasefield_value
  + eta / delta_time * (phasefield_value - old_phasefield)
  + degradation_function_derivative(phasefield_value) * history_value )
  * N_phasefield[i]
  ) * JxW;
  }
  else
  Assert(i_group <= m_d_dof, ExcInternalError());
  }
  }
 

if there is surface pressure, this surface pressure always applied to the reference configuration

  const unsigned int face_pressure_id = 100;
  const double p0 = 0.0;
 
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() && face->boundary_id() == face_pressure_id)
  {
  fe_face_values.reinit(cell, face);
 
  for (const unsigned int f_q_point : fe_face_values.quadrature_point_indices())
  {
  const Tensor<1, dim> &N = fe_face_values.normal_vector(f_q_point);
 
  const double pressure = p0 * time_ramp;
  const Tensor<1, dim> traction = pressure * N;
 
  for (const unsigned int i : fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  const unsigned int component_i = m_fe.system_to_component_index(i).first;
  const double Ni = fe_face_values.shape_value(i, f_q_point);
  const double JxW = fe_face_values.JxW(f_q_point);
  cell_rhs(i) -= (Ni * traction[component_i]) * JxW;
  }
  }
  }
  }
  }
 
  cell->get_dof_indices(local_dof_indices);
  for (const unsigned int i : fe_values.dof_indices())
  system_rhs(local_dof_indices[i]) += cell_rhs(i);
  } // for (const auto &cell : m_dof_handler.active_cell_iterators())
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::update_history_field_step()
  {
  m_logfile << "\t\tUpdate history variable" << std::endl;
 
  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
  {
  lqph[q_point]->update_history_variable();
  }
  }
  }
 
  template <int dim>
  double PhaseFieldMonolithicSolve<dim>::line_search_stepsize_gradient_based(const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta)
  {
  BlockVector<double> g_old(m_system_rhs);
 

BFGS_p_vector is the search direction

  BlockVector<double> solution_delta_trial(solution_delta);

take a full step size 1.0

  solution_delta_trial.add(1.0, BFGS_p_vector);
 
  update_qph_incremental(solution_delta_trial, m_solution, false);
 
  BlockVector<double> g_new(m_dofs_per_block);
  assemble_system_rhs_BFGS_parallel(m_solution, g_new);
 
  BlockVector<double> y_old(m_dofs_per_block);
 
  y_old = g_new - g_old;
 
  double alpha = 1.0;
 
  double alpha_old = 0.0;
 
  double delta_alpha_old = alpha - alpha_old;
 
  double delta_alpha_new;
 
  unsigned int ls_max = 10;
 
  for (unsigned int i = 1; i <= ls_max; ++i)
  {
  delta_alpha_new = -delta_alpha_old
  * (g_new * BFGS_p_vector)/(y_old * BFGS_p_vector);
  alpha += delta_alpha_new;
 
  if (std::fabs(delta_alpha_new) < 1.0e-5)
  break;
 
  if (i == ls_max)
  {
  alpha = 1.0;
  break;
  }
 
  g_old = g_new;
 

BFGS_p_vector is the search direction

  solution_delta_trial = solution_delta;
  solution_delta_trial.add(alpha, BFGS_p_vector);
  update_qph_incremental(solution_delta_trial, m_solution, false);
  assemble_system_rhs_BFGS_parallel(m_solution, g_new);
 
  y_old = g_new - g_old;
 
  delta_alpha_old = delta_alpha_new;
  }
 
  if (alpha < 1.0e-3)
  alpha = 1.0;
 
  return alpha;
  }
 
  template <int dim>
  double PhaseFieldMonolithicSolve<dim>::line_search_stepsize_strong_wolfe(const double phi_0,
  const double phi_0_prime,
  const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta)
  {

AssertThrow(phi_0_prime < 0, ExcMessage("The derivative of phi at alpha = 0 should be negative!"));

Some line search parameters

  const double c1 = 0.0001;
  const double c2 = 0.9;
  const double alpha_max = 100.0;
  const unsigned int max_iter = 20;
  double alpha = 1.0;
 
  double phi_old = phi_0;
  double phi_prime_old = phi_0_prime;
  double alpha_old = 0.0;
 
  double phi, phi_prime;
 
  std::pair<double, double> current_phi_phi_prime;
 
  unsigned int i = 0;
  for (; i < max_iter; ++i)
  {
  current_phi_phi_prime = calculate_phi_and_phi_prime(alpha, BFGS_p_vector, solution_delta);
  phi = current_phi_phi_prime.first;
  phi_prime = current_phi_phi_prime.second;
 
  if ( ( phi > (phi_0 + c1 * alpha * phi_0_prime) )
  || ( i > 0 && phi > phi_old ) )
  {
  return line_search_zoom_strong_wolfe(phi_old, phi_prime_old, alpha_old,
  phi, phi_prime, alpha,
  phi_0, phi_0_prime, BFGS_p_vector,
  c1, c2, max_iter, solution_delta);
  }
 
  if (std::fabs(phi_prime) <= c2 * std::fabs(phi_0_prime))
  {
  return alpha;
  }
 
  if (phi_prime >= 0)
  {
  return line_search_zoom_strong_wolfe(phi, phi_prime, alpha,
  phi_old, phi_prime_old, alpha_old,
  phi_0, phi_0_prime, BFGS_p_vector,
  c1, c2, max_iter, solution_delta);
  }
 
  phi_old = phi;
  phi_prime_old = phi_prime;
  alpha_old = alpha;
 
  alpha = std::min(2.0*alpha, alpha_max);
 
std::min
::VectorizedArray< Number, width > min(const ::VectorizedArray< Number, width > &, const ::VectorizedArray< Number, width > &)
Definition vectorization.h:6596

AssertThrow(alpha < alpha_max, ExcMessage("alpha is bigger than alpha_max, line search failed!"));

  }
 

AssertThrow(i < max_iter, ExcMessage("max number attempts arrived, line search failed!")); Instead of terminating the program, we can just take a full step.

  if (i == max_iter)
  alpha = 1.0;
 
  return alpha;
  }
 
  template <int dim>
  double PhaseFieldMonolithicSolve<dim>::
  line_search_zoom_strong_wolfe(double phi_low, double phi_low_prime, double alpha_low,
  double phi_high, double phi_high_prime, double alpha_high,
  double phi_0, double phi_0_prime, const BlockVector<double> & BFGS_p_vector,
  double c1, double c2, unsigned int max_iter, const BlockVector<double> & solution_delta)
  {
  double alpha = 0;
  std::pair<double, double> current_phi_phi_prime;
  double phi, phi_prime;
 
  unsigned int i = 0;
  for (; i < max_iter; ++i)
  {

a simple bisection is faster than cubic interpolation

  alpha = 0.5 * (alpha_low + alpha_high);

alpha = line_search_interpolation_cubic(alpha_low, phi_low, phi_low_prime, alpha_high, phi_high, phi_high_prime);

  current_phi_phi_prime = calculate_phi_and_phi_prime(alpha, BFGS_p_vector, solution_delta);
  phi = current_phi_phi_prime.first;
  phi_prime = current_phi_phi_prime.second;
 
  if ( (phi > phi_0 + c1 * alpha * phi_0_prime)
  || (phi > phi_low) )
  {
  alpha_high = alpha;
  phi_high = phi;
  phi_high_prime = phi_prime;
  }
  else
  {
  if (std::fabs(phi_prime) <= c2 * std::fabs(phi_0_prime))
  {

if (alpha < 1.0e-3) alpha = 1.0e-3;

  return alpha;
  }
 
  if (phi_prime * (alpha_high - alpha_low) >= 0.0)
  {
  alpha_high = alpha_low;
  phi_high_prime = phi_low_prime;
  phi_high = phi_low;
  }
 
  alpha_low = alpha;
  phi_low_prime = phi_prime;
  phi_low = phi;
  }
  }
 
  if (alpha < 1.0e-3)
  alpha = 1.0;
 

avoid unused variable warnings from compiler

  (void)phi_high;
  (void)phi_high_prime;
  return alpha;
  }
 
  template <int dim>
  double PhaseFieldMonolithicSolve<dim>::
  line_search_interpolation_cubic(const double alpha_0, const double phi_0, const double phi_0_prime,
  const double alpha_1, const double phi_1, const double phi_1_prime)
  {
  const double d1 = phi_0_prime + phi_1_prime - 3.0 * (phi_0 - phi_1) / (alpha_0 - alpha_1);
 
  const double temp = d1 * d1 - phi_0_prime * phi_1_prime;
 
  if (temp < 0.0)
  return 0.5 * (alpha_0 + alpha_1);
 
  int sign;
  if (alpha_1 > alpha_0)
  sign = 1;
  else
  sign = -1;
 
  const double d2 = sign * std::sqrt(temp);
 
  const double alpha = alpha_1 - (alpha_1 - alpha_0)
  * (phi_1_prime + d2 - d1) / (phi_1_prime - phi_0_prime + 2*d2);
 
  if ( (alpha_1 > alpha_0)
  && (alpha > alpha_1 || alpha < alpha_0))
  return 0.5 * (alpha_0 + alpha_1);
 
  if ( (alpha_0 > alpha_1)
  && (alpha > alpha_0 || alpha < alpha_1))
  return 0.5 * (alpha_0 + alpha_1);
 
  return alpha;
  }
 
  template <int dim>
  std::pair<double, double> PhaseFieldMonolithicSolve<dim>::
  calculate_phi_and_phi_prime(const double alpha,
  const BlockVector<double> & BFGS_p_vector,
  const BlockVector<double> & solution_delta)
  {

the first component is phi(alpha), the second component is phi_prime(alpha),

  std::pair<double, double> phi_values;
 
  BlockVector<double> solution_delta_trial(solution_delta);
  solution_delta_trial.add(alpha, BFGS_p_vector);
 
  update_qph_incremental(solution_delta_trial, m_solution, false);
 
  BlockVector<double> system_rhs(m_dofs_per_block);
  assemble_system_rhs_BFGS_parallel(m_solution, system_rhs);

m_constraints.condense(system_rhs);

  phi_values.first = calculate_energy_functional();
  phi_values.second = system_rhs * BFGS_p_vector;
  return phi_values;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::LBFGS_B0(BlockVector<double> & LBFGS_r_vector,
  BlockVector<double> & LBFGS_q_vector)
  {
  m_timer.enter_subsection("Solve B0");
 
  assemble_system_B0(m_solution);
 
  if (m_parameters.m_type_linear_solver == "Direct")
  {
  SparseDirectUMFPACK A_direct;
  A_direct.initialize(m_tangent_matrix);
  A_direct.vmult(LBFGS_r_vector,
  LBFGS_q_vector);
  }
  else if (m_parameters.m_type_linear_solver == "CG")
  {
  /*
  SolverControl solver_control(1e6, 1e-9);
  SolverCG<BlockVector<double>> cg(solver_control);
 
  PreconditionJacobi<BlockSparseMatrix<double>> preconditioner;
  preconditioner.initialize(m_tangent_matrix, 1.0);
 
  cg.solve(m_tangent_matrix,
  LBFGS_r_vector,
  LBFGS_q_vector,
  preconditioner);
  */
  SolverControl solver_control_uu(1e6, 1e-9);
  SolverCG<Vector<double>> cg_uu(solver_control_uu);
 
  PreconditionJacobi<SparseMatrix<double>> preconditioner_uu;
  preconditioner_uu.initialize(m_tangent_matrix.block(m_u_dof, m_u_dof), 1.0);
  cg_uu.solve(m_tangent_matrix.block(m_u_dof, m_u_dof),
  LBFGS_r_vector.block(m_u_dof),
  LBFGS_q_vector.block(m_u_dof),
  preconditioner_uu);
 
  SolverControl solver_control_dd(1e6, 1e-15);
  SolverCG<Vector<double>> cg_dd(solver_control_dd);
 
  PreconditionJacobi<SparseMatrix<double>> preconditioner_dd;
  preconditioner_dd.initialize(m_tangent_matrix.block(m_d_dof, m_d_dof), 1.0);
  cg_dd.solve(m_tangent_matrix.block(m_d_dof, m_d_dof),
  LBFGS_r_vector.block(m_d_dof),
  LBFGS_q_vector.block(m_d_dof),
  preconditioner_dd);
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected linear solver not implemented!"));
  }
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  std::vector<double>
  PhaseFieldMonolithicSolve<dim>::solve_linear_system(BlockVector<double> & newton_update)
  {
  m_timer.enter_subsection("Solve coupled linear system");
 
  if (m_parameters.m_output_iteration_history)
  m_logfile << " SLV " << std::flush;
 
  std::vector<double> linear_solver_parameters(3);
  /*
  {
  SolverControl solver_control(1e6, 1e-9);
  SolverCG<Vector<double>> cg(solver_control);
  cg.connect_condition_number_slot(
  [&] (double condition_number)
  {
  linear_solver_parameters[0] = condition_number;
BlockVectorBase::block
BlockType & block(const unsigned int i)
PreconditionJacobi::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())
SparseDirectUMFPACK::initialize
void initialize(const SparsityPattern &sparsity_pattern)
Definition sparse_direct.cc:74

m_logfile << " Estimated condition number = "<< condition_number << std::endl;

  },
  false);
 
  PreconditionSSOR<SparseMatrix<double>> preconditioner;
  preconditioner.initialize(m_system_matrix_displacement, 1.2);
 
  cg.solve(m_system_matrix_displacement,
  newton_update,
  m_system_rhs_displacement,
  preconditioner);
 
PreconditionSSOR::initialize
void initialize(const MatrixType &A, const AdditionalData &parameters=AdditionalData())

m_logfile << " " << solver_control.last_step() << " CG iterations needed to obtain convergence." << std::endl;

  linear_solver_parameters[1] = solver_control.last_step();
  linear_solver_parameters[2] = solver_control.last_value();
  }
  */
  SparseDirectUMFPACK A_direct;
  A_direct.initialize(m_tangent_matrix);
  A_direct.vmult(newton_update,
  m_system_rhs);
 
  m_constraints.distribute(newton_update);
 
  m_timer.leave_subsection();
  return linear_solver_parameters;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::print_conv_header_newton()
  {
  static const unsigned int l_width = 135;
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
 
  m_logfile << " SOLVER STEP (Newton) "
  << " | Cond No. Lin_Iter Lin_Res Res_Norm "
  << " Res_u Res_d Inc_Norm "
  << " Inc_u Inc_d" << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::print_conv_header_BFGS()
  {
  static const unsigned int l_width = 125;
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
 
  m_logfile << " SOLVER STEP (BFGS) "
  << " | Line Search alpha Energy Res_Norm "
  << " Res_u Res_d Inc_Norm "
  << " Inc_u Inc_d" << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::print_conv_header_LBFGS()
  {
  static const unsigned int l_width = 120;
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
 
  m_logfile << " SOLVER STEP (LBFGS) "
  << " | LS-alpha Energy Res_Norm "
  << " Res_u Res_d Inc_Norm "
  << " Inc_u Inc_d" << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < l_width; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
  }
 
  template <int dim>
  bool PhaseFieldMonolithicSolve<dim>::
  solve_nonlinear_timestep_newton(BlockVector<double> & solution_delta)
  {
  BlockVector<double> newton_update(m_dofs_per_block);
 
  m_error_residual.reset();
  m_error_residual_0.reset();
  m_error_residual_norm.reset();
  m_error_update.reset();
  m_error_update_0.reset();
  m_error_update_norm.reset();
 
  if (m_parameters.m_output_iteration_history)
  print_conv_header_newton();
 
  unsigned int newton_iteration = 0;
  for (; newton_iteration < m_parameters.m_max_iterations_NR; ++newton_iteration)
  {
  if (m_parameters.m_output_iteration_history)
  m_logfile << '\t' << '\t' << std::setw(2) << newton_iteration << ' '
  << std::flush;
 
  make_constraints(newton_iteration);
  assemble_system_newton(m_solution);
 
  get_error_residual(m_error_residual);
  if (newton_iteration == 0)
  m_error_residual_0 = m_error_residual;
 
  m_error_residual_norm = m_error_residual;
  m_error_residual_norm.normalize(m_error_residual_0);
 
  if (newton_iteration > 0 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
  && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
  && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
  && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual)
  {
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " CONVERGED!";
  m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
  << std::scientific
  << " " << " ---- "
  << " " << " ---- "
  << " " << " ---- "
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < 135; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
  }
 
  m_logfile << "\t\tConvergence is reached after "
  << newton_iteration << " Newton iterations."<< std::endl;
 
  m_logfile << "\t\tResidual information of convergence:" << std::endl;
 
  m_logfile << "\t\t\tRelative residual of disp. equation: "
  << m_error_residual_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute residual of disp. equation: "
  << m_error_residual_norm.m_u * m_error_residual_0.m_u << std::endl;
 
  m_logfile << "\t\t\tRelative residual of phasefield equation: "
  << m_error_residual_norm.m_d << std::endl;
 
  m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
  << m_error_residual_norm.m_d * m_error_residual_0.m_d << std::endl;
 
  m_logfile << "\t\t\tRelative increment of disp.: "
  << m_error_update_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of disp.: "
  << m_error_update_norm.m_u * m_error_update_0.m_u << std::endl;
 
  m_logfile << "\t\t\tRelative increment of phasefield: "
  << m_error_update_norm.m_d << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of phasefield: "
  << m_error_update_norm.m_d * m_error_update_0.m_d << std::endl;
 

break;

  return true;
  }
 
  std::vector<double> linear_solver_parameters(3);
 
  linear_solver_parameters = solve_linear_system(newton_update);
 
  get_error_update(newton_update, m_error_update);
  if (newton_iteration == 0)
  m_error_update_0 = m_error_update;
 
  m_error_update_norm = m_error_update;
  m_error_update_norm.normalize(m_error_update_0);
 
  solution_delta += newton_update;
  update_qph_incremental(solution_delta, m_solution, true);
 
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
  << std::scientific
  << " " << linear_solver_parameters[0]
  << " " << linear_solver_parameters[1]
  << " " << linear_solver_parameters[2]
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
  }
  }
 

AssertThrow(newton_iteration < m_parameters.m_max_iterations_NR, ExcMessage("No convergence in Newton-Raphson nonlinear solver!"));

  return false;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::
  solve_nonlinear_timestep_BFGS(BlockVector<double> & solution_delta)
  {
  AssertThrow(false,
  ExcMessage("BFGS requires too much memory. Please use L-BFGS!"));
 
  BlockVector<double> BFGS_update(m_dofs_per_block);
 
  m_error_residual.reset();
  m_error_residual_0.reset();
  m_error_residual_norm.reset();
  m_error_update.reset();
  m_error_update_0.reset();
  m_error_update_norm.reset();
 
  print_conv_header_BFGS();
 
  unsigned int BFGS_iteration = 0;
 

Initial guess B_0, which is a full matrix and takes a lot of memory

  FullMatrix<double> BFGS_matrix = IdentityMatrix(m_dof_handler.n_dofs());
  Vector<double> BFGS_r_vector(m_dof_handler.n_dofs());
  Vector<double> BFGS_p_vector(m_dof_handler.n_dofs());
  Vector<double> BFGS_y_vector(m_dof_handler.n_dofs());
  Vector<double> BFGS_temp_vector(m_dof_handler.n_dofs());
 
  double line_search_parameter, rho;
 

Most likely, we will not be able to create a second full matrix since we will run out of memory on a laptop workstation

  FullMatrix<double> temp_matrix_1(m_dof_handler.n_dofs());
  FullMatrix<double> temp_matrix_2(m_dof_handler.n_dofs());
 
  for (; BFGS_iteration < m_parameters.m_max_iterations_BFGS; ++BFGS_iteration)
  {
  m_logfile << '\t' << '\t' << std::setw(2) << BFGS_iteration << ' '
  << std::flush;
 
  make_constraints(BFGS_iteration);
 

At the first step, we simply distribute the inhomogeneous part of the constraints

  if (BFGS_iteration == 0)
  {
  m_constraints.distribute(BFGS_update);
  solution_delta += BFGS_update;
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
  update_qph_incremental(solution_delta, m_solution, false);
  m_logfile << " --- |" << std::flush;
  m_logfile << std::endl;
  continue;
  }
  else if (BFGS_iteration == 1)
  {

Calculate the residual vector r. NOTICE that in the context of BFGS, this r is the gradient of the energy functional (objective function), NOT the negative gradient of the energy functional

  assemble_system_rhs_BFGS(m_solution, m_system_rhs);
 

We cannot simply zero out the dofs that are constrained, since we might have hanging node constraints. In this case, we need to modify the RHS as C^T * b, which C contains entries of 0.5 (x_3 = 0.5*x_1 + 0.5*x_2) for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i) if (m_constraints.is_constrained(i)) m_system_rhs(i) = 0.0;

if m_constraints has inhomogeneity, we cannot call m_constraints.condense(m_system_rhs), since the m_system_matrix needs to be provided to modify the RHS properly. However, this error will not be detected in the release mode and only will be detected on the debug mode

  m_constraints.condense(m_system_rhs);
  }
 
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
 
  get_error_residual(m_error_residual);
  if (BFGS_iteration == 1)
  m_error_residual_0 = m_error_residual;
 
  m_error_residual_norm = m_error_residual;
  m_error_residual_norm.normalize(m_error_residual_0);
 
  if (BFGS_iteration > 1 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
  && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
  && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
  && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual)
  {
  m_logfile << " CONVERGED!";
  m_logfile << "| " << std::fixed << std::setprecision(3) << std::setw(7)
  << std::scientific
  << " " << " ---- "
  << " " << " ---- "
  << " " << " ---- "
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < 135; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
 
  break;
  }
 

BFGS algorithm

  BFGS_r_vector = m_system_rhs;
  BFGS_matrix.vmult(BFGS_p_vector, BFGS_r_vector);
  BFGS_p_vector *= -1.0;
  m_constraints.distribute(BFGS_p_vector);
 

We need a line search algorithm to decide line_search_parameter

  const double phi_0 = calculate_energy_functional();
  const double phi_0_prime = BFGS_r_vector * BFGS_p_vector;
 
  BlockVector<double> BFGS_p_vector_block(m_dofs_per_block);
  BFGS_p_vector_block = BFGS_p_vector;
  line_search_parameter = line_search_stepsize_strong_wolfe(phi_0,
  phi_0_prime,
  BFGS_p_vector_block,
  solution_delta);
 
  BFGS_p_vector *= line_search_parameter;
  BFGS_update = BFGS_p_vector;
 
  get_error_update(BFGS_update, m_error_update);
  if (BFGS_iteration == 1)
  m_error_update_0 = m_error_update;
 
  m_error_update_norm = m_error_update;
  m_error_update_norm.normalize(m_error_update_0);
 
  solution_delta += BFGS_update;
  update_qph_incremental(solution_delta, m_solution, false);
 
  BFGS_y_vector = m_system_rhs;
  BFGS_y_vector *= -1.0;
  assemble_system_rhs_BFGS(m_solution, m_system_rhs);
  m_constraints.condense(m_system_rhs);
  BFGS_temp_vector = m_system_rhs;
  BFGS_y_vector += BFGS_temp_vector;
 

rho should be positive with the proper line search

  rho = BFGS_y_vector * BFGS_p_vector;
  rho = 1.0/rho;
 
  if (rho < 0)
  m_logfile << "Rho is negative!" << std::endl;
 

In the first step, we scale the identity matrix as the BFGS matrix

  if (BFGS_iteration == 1)
  {
  double scale_parameter = (BFGS_y_vector * BFGS_p_vector) / (BFGS_y_vector.norm_sqr());
  BFGS_matrix *= scale_parameter;
  }
 
  temp_matrix_1.outer_product(BFGS_p_vector, BFGS_y_vector);
  temp_matrix_2 = IdentityMatrix(m_dof_handler.n_dofs());
  temp_matrix_2.add(-rho, temp_matrix_1);
 
  temp_matrix_2.mmult(temp_matrix_1, BFGS_matrix);
  temp_matrix_1.mTmult(BFGS_matrix, temp_matrix_2);
 
  temp_matrix_1.outer_product(BFGS_p_vector, BFGS_p_vector);
 
  BFGS_matrix.add(rho, temp_matrix_1);
 
  const double energy_functional = calculate_energy_functional();
 
  m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(7)
  << std::scientific
  << " " << line_search_parameter
  << " " << energy_functional
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
  }
 
  AssertThrow(BFGS_iteration < m_parameters.m_max_iterations_BFGS,
  ExcMessage("No convergence in BFGS nonlinear solver!"));
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::
  solve_nonlinear_timestep_LBFGS(BlockVector<double> & solution_delta,
  BlockVector<double> & LBFGS_update_refine)
  {
  BlockVector<double> LBFGS_update(m_dofs_per_block);
 
  LBFGS_update = 0.0;
 
  m_error_residual.reset();
  m_error_residual_0.reset();
  m_error_residual_norm.reset();
  m_error_update.reset();
  m_error_update_0.reset();
  m_error_update_norm.reset();
 
  if (m_parameters.m_output_iteration_history)
  print_conv_header_LBFGS();
 
  unsigned int LBFGS_iteration = 0;
 
  BlockVector<double> LBFGS_r_vector(m_dofs_per_block);
  BlockVector<double> LBFGS_y_vector(m_dofs_per_block);
  BlockVector<double> LBFGS_q_vector(m_dofs_per_block);
  BlockVector<double> LBFGS_s_vector(m_dofs_per_block);
  std::list<std::pair< std::pair<BlockVector<double>,
  BlockVector<double>>,
  double>> LBFGS_vector_list;
 
  const unsigned int LBFGS_m = m_parameters.m_LBFGS_m;
  std::list<double> LBFGS_alpha_list;
 
  double line_search_parameter = 0.0;
  double LBFGS_beta = 0.0;
  double rho = 0.0;
 
  for (; LBFGS_iteration < m_parameters.m_max_iterations_BFGS; ++LBFGS_iteration)
  {
  if (m_parameters.m_output_iteration_history)
  m_logfile << '\t' << '\t' << std::setw(2) << LBFGS_iteration << ' '
  << std::flush;
 
  make_constraints(LBFGS_iteration);
 

At the first step, we simply distribute the inhomogeneous part of the constraints

  if (LBFGS_iteration == 0)
  {

use the solution from the previous solve on the refined mesh as initial guess

  LBFGS_update = LBFGS_update_refine;
 
  m_constraints.distribute(LBFGS_update);
  solution_delta += LBFGS_update;
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
  }
  update_qph_incremental(solution_delta, m_solution, false);
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " --- |" << std::flush;
  m_logfile << std::endl;
  }
  continue;
  }
  else if (LBFGS_iteration == 1)
  {

Calculate the residual vector r. NOTICE that in the context of BFGS, this r is the gradient of the energy functional (objective function), NOT the negative gradient of the energy functional

  assemble_system_rhs_BFGS_parallel(m_solution, m_system_rhs);
 

We cannot simply zero out the dofs that are constrained, since we might have hanging node constraints. In this case, we need to modify the RHS as C^T * b, which C contains entries of 0.5 (x_3 = 0.5*x_1 + 0.5*x_2) for (unsigned int i = 0; i < m_dof_handler.n_dofs(); ++i) if (m_constraints.is_constrained(i)) m_system_rhs(i) = 0.0;

if m_constraints has inhomogeneity, we cannot call m_constraints.condense(m_system_rhs), since the m_system_matrix needs to be provided to modify the RHS properly. However, this error will not be detected in the release mode and only will be detected on the debug mode if we use assemble_system_rhs_BFGS_parallel, then condense() is not necessary m_constraints.condense(m_system_rhs);

  }
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
  m_logfile << " --- " << std::flush;
  }
 
  get_error_residual(m_error_residual);
  if (LBFGS_iteration == 1)
  m_error_residual_0 = m_error_residual;
 
  m_error_residual_norm = m_error_residual;

For three-point bending problem and 3D problem, we use absolute residual for convergence test

  if (m_parameters.m_relative_residual)
  m_error_residual_norm.normalize(m_error_residual_0);
 
  if (LBFGS_iteration > 1 && m_error_update_norm.m_u <= m_parameters.m_tol_u_incr
  && m_error_residual_norm.m_u <= m_parameters.m_tol_u_residual
  && m_error_update_norm.m_d <= m_parameters.m_tol_d_incr
  && m_error_residual_norm.m_d <= m_parameters.m_tol_d_residual
  )
  {
  if (m_parameters.m_output_iteration_history)
  {
  m_logfile << " | ";
  m_logfile << " CONVERGED! " << std::fixed << std::setprecision(3) << std::setw(7)
  << std::scientific
  << " ---- "
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
 
  m_logfile << '\t' << '\t';
  for (unsigned int i = 0; i < 120; ++i)
  m_logfile << '_';
  m_logfile << std::endl;
  }
 
  m_logfile << "\t\tConvergence is reached after "
  << LBFGS_iteration << " L-BFGS iterations."<< std::endl;
 
  m_logfile << "\t\tResidual information of convergence:" << std::endl;
 
  if (m_parameters.m_relative_residual)
  {
  m_logfile << "\t\t\tRelative residual of disp. equation: "
  << m_error_residual_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute residual of disp. equation: "
  << m_error_residual_norm.m_u * m_error_residual_0.m_u << std::endl;
 
  m_logfile << "\t\t\tRelative residual of phasefield equation: "
  << m_error_residual_norm.m_d << std::endl;
 
  m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
  << m_error_residual_norm.m_d * m_error_residual_0.m_d << std::endl;
 
  m_logfile << "\t\t\tRelative increment of disp.: "
  << m_error_update_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of disp.: "
  << m_error_update_norm.m_u * m_error_update_0.m_u << std::endl;
 
  m_logfile << "\t\t\tRelative increment of phasefield: "
  << m_error_update_norm.m_d << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of phasefield: "
  << m_error_update_norm.m_d * m_error_update_0.m_d << std::endl;
  }
  else
  {
  m_logfile << "\t\t\tAbsolute residual of disp. equation: "
  << m_error_residual_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute residual of phasefield equation: "
  << m_error_residual_norm.m_d << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of disp.: "
  << m_error_update_norm.m_u << std::endl;
 
  m_logfile << "\t\t\tAbsolute increment of phasefield: "
  << m_error_update_norm.m_d << std::endl;
  }
 
  break;
  }
 

LBFGS algorithm

  LBFGS_q_vector = m_system_rhs;
 
  LBFGS_alpha_list.clear();
  for (auto itr = LBFGS_vector_list.begin(); itr != LBFGS_vector_list.end(); ++itr)
  {
  LBFGS_s_vector = (itr->first).first;
  LBFGS_y_vector = (itr->first).second;
  rho = itr->second;
 
  const double alpha = rho * (LBFGS_s_vector * LBFGS_q_vector);
  LBFGS_alpha_list.push_back(alpha);
 
  LBFGS_q_vector.add(-alpha, LBFGS_y_vector);
  }
  /*
  double scale_gamma = 0.0;
  if (LBFGS_iteration == 1)
  {
  scale_gamma = 1.0;
  }
  else
  {
  LBFGS_s_vector = LBFGS_vector_list.front().first.first;
  LBFGS_y_vector = LBFGS_vector_list.front().first.second;
  scale_gamma = (LBFGS_s_vector * LBFGS_y_vector)/(LBFGS_y_vector * LBFGS_y_vector);
  }
 
  LBFGS_q_vector *= scale_gamma;
  LBFGS_r_vector = LBFGS_q_vector;
  */
  LBFGS_B0(LBFGS_r_vector,
  LBFGS_q_vector);
 
  for (auto itr = LBFGS_vector_list.rbegin(); itr != LBFGS_vector_list.rend(); ++itr)
  {
  LBFGS_s_vector = (itr->first).first;
  LBFGS_y_vector = (itr->first).second;
  rho = itr->second;
 
  LBFGS_beta = rho * (LBFGS_y_vector * LBFGS_r_vector);
 
  const double alpha = LBFGS_alpha_list.back();
  LBFGS_alpha_list.pop_back();
 
  LBFGS_r_vector.add(alpha - LBFGS_beta, LBFGS_s_vector);
  }
 
  LBFGS_r_vector *= -1.0; // this is the p_vector (search direction)
 
  m_constraints.distribute(LBFGS_r_vector);
 
BlockVectorBase::add
void add(const std::vector< size_type > &indices, const std::vector< Number > &values)
second
Point< 2 > second
Definition grid_out.cc:4640

We need a line search algorithm to decide line_search_parameter

  if(m_parameters.m_type_line_search == "StrongWolfe")
  {
  const double phi_0 = calculate_energy_functional();
  const double phi_0_prime = m_system_rhs * LBFGS_r_vector;
 
  line_search_parameter = line_search_stepsize_strong_wolfe(phi_0,
  phi_0_prime,
  LBFGS_r_vector,
  solution_delta);
  }
  else if(m_parameters.m_type_line_search == "GradientBased")
  {

LBFGS_r_vector is the search direction

  line_search_parameter = line_search_stepsize_gradient_based(LBFGS_r_vector,
  solution_delta);
  }
  else
  {
  Assert(false, ExcMessage("An unknown line search method is called!"));
  }
 
  LBFGS_r_vector *= line_search_parameter;
  LBFGS_update = LBFGS_r_vector;
 
  get_error_update(LBFGS_update, m_error_update);
  if (LBFGS_iteration == 1)
  m_error_update_0 = m_error_update;
 
  m_error_update_norm = m_error_update;

For three-point bending problem and the sphere inclusion problem, we use absolute residual for convergence test

  if (m_parameters.m_relative_residual)
  m_error_update_norm.normalize(m_error_update_0);
 
  solution_delta += LBFGS_update;
  update_qph_incremental(solution_delta, m_solution, false);
 
  LBFGS_y_vector = m_system_rhs;
  LBFGS_y_vector *= -1.0;
  assemble_system_rhs_BFGS_parallel(m_solution, m_system_rhs);

if we use assemble_system_rhs_BFGS_parallel, then condense() is not necessary m_constraints.condense(m_system_rhs);

  LBFGS_y_vector += m_system_rhs;
 
  LBFGS_s_vector = LBFGS_update;
 
  const double g_norm = m_system_rhs.l2_norm();
 
  const double yxs = LBFGS_y_vector * LBFGS_s_vector;
 
  const double sxs = LBFGS_s_vector * LBFGS_s_vector;
 
  if (yxs/sxs >= 1.0e-6 * g_norm)
  {
  if (LBFGS_iteration > LBFGS_m)
  LBFGS_vector_list.pop_back();
 
  rho = 1.0 / yxs;
 
  LBFGS_vector_list.push_front(std::make_pair(std::make_pair(LBFGS_s_vector,
  LBFGS_y_vector),
  rho));
  }
 
  if (m_parameters.m_output_iteration_history)
  {
  const double energy_functional = calculate_energy_functional();
 
  m_logfile << " | " << std::fixed << std::setprecision(3) << std::setw(1)
  << std::scientific
  << "" << line_search_parameter
  << std::fixed << std::setprecision(6) << std::setw(1)
  << std::scientific
  << " " << energy_functional
  << std::fixed << std::setprecision(3) << std::setw(1)
  << std::scientific
  << " " << m_error_residual_norm.m_norm
  << " " << m_error_residual_norm.m_u
  << " " << m_error_residual_norm.m_d
  << " " << m_error_update_norm.m_norm
  << " " << m_error_update_norm.m_u
  << " " << m_error_update_norm.m_d
  << " " << std::endl;
  }
  }
 
  AssertThrow(LBFGS_iteration < m_parameters.m_max_iterations_BFGS,
  ExcMessage("No convergence in L-BFGS nonlinear solver!"));
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::output_results() const
  {
  m_timer.enter_subsection("Output results");
 
  DataOut<dim> data_out;
 
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
  data_component_interpretation(
  dim, DataComponentInterpretation::component_is_part_of_vector);
 
  data_component_interpretation.push_back(
  DataComponentInterpretation::component_is_scalar);
 
  std::vector<std::string> solution_name(dim, "displacement");
  solution_name.emplace_back("phasefield");
 
  data_out.attach_dof_handler(m_dof_handler);
  data_out.add_data_vector(m_solution,
  solution_name,
  DataOut<dim>::type_dof_data,
  data_component_interpretation);
 
 
  Vector<double> cell_material_id(m_triangulation.n_active_cells());

output material ID for each cell

  for (const auto &cell : m_triangulation.active_cell_iterators())
  {
  cell_material_id(cell->active_cell_index()) = cell->material_id();
  }
  data_out.add_data_vector(cell_material_id, "materialID");
 

Stress L2 projection

  DoFHandler<dim> stresses_dof_handler_L2(m_triangulation);
  FE_Q<dim> stresses_fe_L2(m_parameters.m_poly_degree); //FE_Q element is continuous
  stresses_dof_handler_L2.distribute_dofs(stresses_fe_L2);
  AffineConstraints<double> constraints;
  constraints.clear();
  DoFTools::make_hanging_node_constraints(stresses_dof_handler_L2, constraints);
  constraints.close();
  std::vector<DataComponentInterpretation::DataComponentInterpretation>
  data_component_interpretation_stress(1,
  DataComponentInterpretation::component_is_scalar);
 
  for (unsigned int i = 0; i < dim; ++i)
  for (unsigned int j = i; j < dim; ++j)
  {
  Vector<double> stress_field_L2;
  stress_field_L2.reinit(stresses_dof_handler_L2.n_dofs());
 
  MappingQ<dim> mapping(m_parameters.m_poly_degree + 1);
  VectorTools::project(mapping,
  stresses_dof_handler_L2,
  constraints,
  m_qf_cell,
  [&] (const typename DoFHandler<dim>::active_cell_iterator & cell,
  const unsigned int q) -> double
  {
  return m_quadrature_point_history.get_data(cell)[q]->get_cauchy_stress()[i][j];
  },
  stress_field_L2);
 
  std::string stress_name = "Cauchy_stress_" + std::to_string(i+1) + std::to_string(j+1)
  + "_L2";
 
  data_out.add_data_vector(stresses_dof_handler_L2,
  stress_field_L2,
  stress_name,
  data_component_interpretation_stress);
  }
 
  data_out.build_patches(m_parameters.m_poly_degree);
 
  std::ofstream output("Solution-" + std::to_string(dim) + "d-" +
  Utilities::int_to_string(m_time.get_timestep(),4) + ".vtu");
 
  data_out.write_vtu(output);
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::calculate_reaction_force(unsigned int face_ID)
  {
  m_timer.enter_subsection("Calculate reaction force");
 
  BlockVector<double> system_rhs;
  system_rhs.reinit(m_dofs_per_block);
 
  Vector<double> cell_rhs(m_dofs_per_cell);
  std::vector<types::global_dof_index> local_dof_indices(m_dofs_per_cell);
 
  const double time_ramp = (m_time.current() / m_time.end());
  std::vector<Tensor<1, dim>> rhs_values(m_n_q_points);
  const UpdateFlags uf_cell(update_values | update_gradients |
  update_quadrature_points | update_JxW_values);
  const UpdateFlags uf_face(update_values | update_normal_vectors |
  update_JxW_values);
 
  FEValues<dim> fe_values(m_fe, m_qf_cell, uf_cell);
  FEFaceValues<dim> fe_face_values(m_fe, m_qf_face, uf_face);
 
BlockVector::reinit
void reinit(const unsigned int n_blocks, const size_type block_size=0, const bool omit_zeroing_entries=false)
Vector::reinit
virtual void reinit(const size_type N, const bool omit_zeroing_entries=false)
VectorTools::project
void project(const Mapping< dim, spacedim > &mapping, const DoFHandler< dim, spacedim > &dof, const AffineConstraints< typename VectorType::value_type > &constraints, const Quadrature< dim > &quadrature, const Function< spacedim, typename VectorType::value_type > &function, VectorType &vec, const bool enforce_zero_boundary=false, const Quadrature< dim - 1 > &q_boundary=(dim > 1 ? QGauss< dim - 1 >(2) :Quadrature< dim - 1 >()), const bool project_to_boundary_first=false)

shape function values for displacement field

  std::vector<std::vector<Tensor<1, dim>>>
  Nx(m_qf_cell.size(), std::vector<Tensor<1, dim>>(m_dofs_per_cell));
  std::vector<std::vector<Tensor<2, dim>>>
  grad_Nx(m_qf_cell.size(), std::vector<Tensor<2, dim>>(m_dofs_per_cell));
  std::vector<std::vector<SymmetricTensor<2, dim>>>
  symm_grad_Nx(m_qf_cell.size(), std::vector<SymmetricTensor<2, dim>>(m_dofs_per_cell));
 
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {

if calculate_reaction_force() is defined as const, then we also need to put a const in std::shared_ptr, that is, std::shared_ptr<const PointHistory<dim>>

  const std::vector<std::shared_ptr< PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
  cell_rhs = 0.0;
  fe_values.reinit(cell);
  right_hand_side(fe_values.get_quadrature_points(),
  rhs_values,
  m_parameters.m_x_component*time_ramp,
  m_parameters.m_y_component*time_ramp,
  m_parameters.m_z_component*time_ramp);
 
  for (const unsigned int q_point : fe_values.quadrature_point_indices())
  {
  for (const unsigned int k : fe_values.dof_indices())
  {
  const unsigned int k_group = m_fe.system_to_base_index(k).first.first;
 
  if (k_group == m_u_dof)
  {
  Nx[q_point][k] = fe_values[m_u_fe].value(k, q_point);
  grad_Nx[q_point][k] = fe_values[m_u_fe].gradient(k, q_point);
  symm_grad_Nx[q_point][k] = symmetrize(grad_Nx[q_point][k]);
  }
  }
  }
 
  for (const unsigned int q_point : fe_values.quadrature_point_indices())
  {
  const SymmetricTensor<2, dim> & cauchy_stress = lqph[q_point]->get_cauchy_stress();
 
  const std::vector<Tensor<1,dim>> & N = Nx[q_point];
  const std::vector<SymmetricTensor<2, dim>> & symm_grad_N = symm_grad_Nx[q_point];
  const double JxW = fe_values.JxW(q_point);
 
  for (const unsigned int i : fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  cell_rhs(i) -= (symm_grad_N[i] * cauchy_stress) * JxW;

contributions from the body force to right-hand side

  cell_rhs(i) += N[i] * rhs_values[q_point] * JxW;
  }
  }
  }
 

if there is surface pressure, this surface pressure always applied to the reference configuration

  const unsigned int face_pressure_id = 100;
  const double p0 = 0.0;
 
  for (const auto &face : cell->face_iterators())
  {
  if (face->at_boundary() && face->boundary_id() == face_pressure_id)
  {
  fe_face_values.reinit(cell, face);
 
  for (const unsigned int f_q_point : fe_face_values.quadrature_point_indices())
  {
  const Tensor<1, dim> &N = fe_face_values.normal_vector(f_q_point);
 
  const double pressure = p0 * time_ramp;
  const Tensor<1, dim> traction = pressure * N;
 
  for (const unsigned int i : fe_values.dof_indices())
  {
  const unsigned int i_group = m_fe.system_to_base_index(i).first.first;
 
  if (i_group == m_u_dof)
  {
  const unsigned int component_i = m_fe.system_to_component_index(i).first;
  const double Ni = fe_face_values.shape_value(i, f_q_point);
  const double JxW = fe_face_values.JxW(f_q_point);
  cell_rhs(i) += (Ni * traction[component_i]) * JxW;
  }
  }
  }
  }
  }
 
  cell->get_dof_indices(local_dof_indices);
  for (const unsigned int i : fe_values.dof_indices())
  system_rhs(local_dof_indices[i]) += cell_rhs(i);
  } // for (const auto &cell : m_dof_handler.active_cell_iterators())
 

The difference between the above assembled system_rhs and m_system_rhs is that m_system_rhs is condensed by the m_constraints, which zero out the rhs values associated with the constrained DOFs and modify the rhs values associated with the unconstrained DOFs.

  std::vector< types::global_dof_index > mapping;
  std::set<types::boundary_id> boundary_ids;
  boundary_ids.insert(face_ID);
  DoFTools::map_dof_to_boundary_indices(m_dof_handler,
  boundary_ids,
  mapping);
 
  std::vector<double> reaction_force(dim, 0.0);
 
  for (unsigned int i = 0; i < m_dofs_per_block[m_u_dof]; ++i)
  {
  if (mapping[i] != numbers::invalid_dof_index)
  {
  reaction_force[i % dim] += system_rhs.block(m_u_dof)(i);
  }
  }
 
  for (unsigned int i = 0; i < dim; i++)
  m_logfile << "\t\tReaction force in direction " << i << " on boundary ID " << face_ID
  << " = "
  << std::fixed << std::setprecision(3) << std::setw(1)
  << std::scientific
  << reaction_force[i] << std::endl;
 
  std::pair<double, std::vector<double>> time_force;
  time_force.first = m_time.current();
  time_force.second = reaction_force;
  m_history_reaction_force.push_back(time_force);
 
  m_timer.leave_subsection();
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::write_history_data()
  {
  m_logfile << "\t\tWrite history data ... \n"<<std::endl;
 
  std::ofstream myfile_reaction_force ("Reaction_force.hist");
  if (myfile_reaction_force.is_open())
  {
  myfile_reaction_force << 0.0 << "\t";
  if (dim == 2)
  myfile_reaction_force << 0.0 << "\t"
  << 0.0 << std::endl;
  if (dim == 3)
  myfile_reaction_force << 0.0 << "\t"
  << 0.0 << "\t"
  << 0.0 << std::endl;
 
  for (const auto &time_force : m_history_reaction_force)
  {
  myfile_reaction_force << time_force.first << "\t";
  if (dim == 2)
  myfile_reaction_force << time_force.second[0] << "\t"
  << time_force.second[1] << std::endl;
  if (dim == 3)
  myfile_reaction_force << time_force.second[0] << "\t"
  << time_force.second[1] << "\t"
  << time_force.second[2] << std::endl;
  }
  myfile_reaction_force.close();
  }
  else
  m_logfile << "Unable to open file";
 
  std::ofstream myfile_energy ("Energy.hist");
  if (myfile_energy.is_open())
  {
  myfile_energy << std::fixed << std::setprecision(10) << std::scientific
  << 0.0 << "\t"
  << 0.0 << "\t"
  << 0.0 << "\t"
  << 0.0 << std::endl;
 
  for (const auto &time_energy : m_history_energy)
  {
  myfile_energy << std::fixed << std::setprecision(10) << std::scientific
  << time_energy.first << "\t"
  << time_energy.second[0] << "\t"
  << time_energy.second[1] << "\t"
  << time_energy.second[2] << std::endl;
  }
  myfile_energy.close();
  }
  else
  m_logfile << "Unable to open file";
  }
 
  template <int dim>
  double PhaseFieldMonolithicSolve<dim>::calculate_energy_functional() const
  {
  double energy_functional = 0.0;
 
  FEValues<dim> fe_values(m_fe, m_qf_cell, update_JxW_values);
 
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
 
  const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
  {
  const double JxW = fe_values.JxW(q_point);
  energy_functional += lqph[q_point]->get_total_strain_energy() * JxW;
  energy_functional += lqph[q_point]->get_crack_energy_dissipation() * JxW;
  }
  }
 
  return energy_functional;
  }
 
  template <int dim>
  std::pair<double, double>
  PhaseFieldMonolithicSolve<dim>::calculate_total_strain_energy_and_crack_energy_dissipation() const
  {
  double total_strain_energy = 0.0;
  double crack_energy_dissipation = 0.0;
 
  FEValues<dim> fe_values(m_fe, m_qf_cell, update_JxW_values);
 
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {
  fe_values.reinit(cell);
 
  const std::vector<std::shared_ptr<const PointHistory<dim>>> lqph =
  m_quadrature_point_history.get_data(cell);
  Assert(lqph.size() == m_n_q_points, ExcInternalError());
 
  for (unsigned int q_point = 0; q_point < m_n_q_points; ++q_point)
  {
  const double JxW = fe_values.JxW(q_point);
  total_strain_energy += lqph[q_point]->get_total_strain_energy() * JxW;
  crack_energy_dissipation += lqph[q_point]->get_crack_energy_dissipation() * JxW;
  }
  }
 
  return std::make_pair(total_strain_energy, crack_energy_dissipation);
  }
 
 
  template <int dim>
  bool PhaseFieldMonolithicSolve<dim>::local_refine_and_solution_transfer(BlockVector<double> & solution_delta,
  BlockVector<double> & LBFGS_update_refine)
  {
push_back
***mech_lbc_system increment_interpolation_handlers push_back(scale_z_handler)
DoFTools::map_dof_to_boundary_indices
void map_dof_to_boundary_indices(const DoFHandler< dim, spacedim > &dof_handler, std::vector< types::global_dof_index > &mapping)
Definition dof_tools.cc:2456
internal::reinit
void reinit(MatrixBlock< MatrixType > &v, const BlockSparsityPattern &p)
Definition matrix_block.h:595
numbers::invalid_dof_index
constexpr types::global_dof_index invalid_dof_index
Definition types.h:259

This is the solution at (n+1) obtained from the old (coarse) mesh

  BlockVector<double> solution_next_step(m_dofs_per_block);
  solution_next_step = m_solution + solution_delta;
  bool mesh_is_same = true;
  bool cell_refine_flag = true;
 
  unsigned int material_id;
  double length_scale;
  double cell_length;
  while(cell_refine_flag)
  {
  cell_refine_flag = false;
 
  std::vector<types::global_dof_index> local_dof_indices(m_fe.dofs_per_cell);
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {
  cell->get_dof_indices(local_dof_indices);
 
  for (unsigned int i = 0; i< m_fe.dofs_per_cell; ++i)
  {
  const unsigned int comp_i = m_fe.system_to_component_index(i).first;
  if (comp_i == m_d_component) //phasefield component
  {
  if ( solution_next_step(local_dof_indices[i])
  > m_parameters.m_phasefield_refine_threshold )
  {
  material_id = cell->material_id();
  length_scale = m_material_data[material_id][2];
  if (dim == 2)
  cell_length = std::sqrt(cell->measure());
  else
  cell_length = std::cbrt(cell->measure());
  if ( cell_length
  > length_scale * m_parameters.m_allowed_max_h_l_ratio )
  {
  if (cell->level() < m_parameters.m_max_allowed_refinement_level)
  {
  cell->set_refine_flag();
  break;
  }
  }
  }
  }
  }
  }
 
  for (const auto &cell : m_dof_handler.active_cell_iterators())
  {
  if (cell->refine_flag_set())
  {
  cell_refine_flag = true;
  break;
  }
  }
 

if any cell is refined, we need to project the solution to the newly refined mesh

  if (cell_refine_flag)
  {
  mesh_is_same = false;
 
  std::vector<BlockVector<double> > old_solutions(2);
  old_solutions[0] = solution_next_step;
  old_solutions[1] = m_solution;
 
  m_triangulation.prepare_coarsening_and_refinement();
  SolutionTransfer<dim, BlockVector<double>> solution_transfer(m_dof_handler);
  solution_transfer.prepare_for_coarsening_and_refinement(old_solutions);
  m_triangulation.execute_coarsening_and_refinement();
 
  setup_system();
 
  std::vector<BlockVector<double>> tmp_solutions(2);
  tmp_solutions[0].reinit(m_dofs_per_block);
  tmp_solutions[1].reinit(m_dofs_per_block);
  solution_transfer.interpolate(tmp_solutions);
  solution_next_step = tmp_solutions[0];
  m_solution = tmp_solutions[1];
 

make sure the projected solutions still satisfy hanging node constraints

  m_constraints.distribute(solution_next_step);
  m_constraints.distribute(m_solution);
  } // if (cell_refine_flag)
  } // while(cell_refine_flag)
 

calculate field variables for newly refined cells

  if (!mesh_is_same)
  {
  BlockVector<double> temp_solution_delta(m_dofs_per_block);
  BlockVector<double> temp_previous_solution(m_dofs_per_block);
  temp_solution_delta = 0.0;
  temp_previous_solution = 0.0;
  update_qph_incremental(temp_solution_delta, temp_previous_solution, false);
  update_history_field_step();
 

initial guess for the resolve on the refined mesh

  LBFGS_update_refine = solution_next_step - m_solution;
  }
 
  return mesh_is_same;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::print_parameter_information()
  {
  m_logfile << "Scenario number = " << m_parameters.m_scenario << std::endl;
  m_logfile << "Log file = " << m_parameters.m_logfile_name << std::endl;
  m_logfile << "Write iteration history to log file? = " << std::boolalpha
  << m_parameters.m_output_iteration_history << std::endl;
  m_logfile << "Nonlinear solver type = " << m_parameters.m_type_nonlinear_solver << std::endl;
  m_logfile << "Line search type = " << m_parameters.m_type_line_search << std::endl;
  m_logfile << "Linear solver type = " << m_parameters.m_type_linear_solver << std::endl;
  m_logfile << "Mesh refinement strategy = " << m_parameters.m_refinement_strategy << std::endl;
  m_logfile << "L-BFGS_m = " << m_parameters.m_LBFGS_m << std::endl;
  m_logfile << "Global refinement times = " << m_parameters.m_global_refine_times << std::endl;
  m_logfile << "Local prerefinement times = " <<m_parameters. m_local_prerefine_times << std::endl;
  m_logfile << "Maximum adaptive refinement times allowed in each step = "
  << m_parameters.m_max_adaptive_refine_times << std::endl;
  m_logfile << "Maximum allowed cell refinement level = "
  << m_parameters.m_max_allowed_refinement_level << std::endl;
  m_logfile << "Phasefield-based refinement threshold value = "
  << m_parameters.m_phasefield_refine_threshold << std::endl;
  m_logfile << "Allowed maximum h/l ratio = " << m_parameters.m_allowed_max_h_l_ratio << std::endl;
  m_logfile << "total number of material types = " << m_parameters.m_total_material_regions << std::endl;
  m_logfile << "material data file name = " << m_parameters.m_material_file_name << std::endl;
  if (m_parameters.m_reaction_force_face_id >= 0)
  m_logfile << "Calculate reaction forces on Face ID = " << m_parameters.m_reaction_force_face_id << std::endl;
  else
  m_logfile << "No need to calculate reaction forces." << std::endl;
 
  if (m_parameters.m_relative_residual)
  m_logfile << "Relative residual for convergence." << std::endl;
  else
  m_logfile << "Absolute residual for convergence." << std::endl;
 
  m_logfile << "Body force = (" << m_parameters.m_x_component << ", "
  << m_parameters.m_y_component << ", "
  << m_parameters.m_z_component << ") (N/m^3)"
  << std::endl;
 
  m_logfile << "End time = " << m_parameters.m_end_time << std::endl;
  m_logfile << "Time data file name = " << m_parameters.m_time_file_name << std::endl;
  }
 
  template <int dim>
  void PhaseFieldMonolithicSolve<dim>::run()
  {
  print_parameter_information();
 
  read_material_data(m_parameters.m_material_file_name,
  m_parameters.m_total_material_regions);
 
  std::vector<std::array<double, 4>> time_table;
 
  read_time_data(m_parameters.m_time_file_name, time_table);
 
  make_grid();
  setup_system();
  output_results();
 
  m_time.increment(time_table);
 
  while(m_time.current() < m_time.end() + m_time.get_delta_t()*1.0e-6)
  {
  m_logfile << std::endl
  << "Timestep " << m_time.get_timestep() << " @ " << m_time.current()
  << 's' << std::endl;
 
  bool mesh_is_same = false;
 

initial guess for the resolve on the refined mesh

  BlockVector<double> LBFGS_update_refine(m_dofs_per_block);
  LBFGS_update_refine = 0.0;
 

local adaptive mesh refinement loop

  unsigned int adp_refine_iteration = 0;
  for (; adp_refine_iteration < m_parameters.m_max_adaptive_refine_times + 1; ++adp_refine_iteration)
  {
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  m_logfile << "\tAdaptive refinement-"<< adp_refine_iteration << ": " << std::endl;
 
  BlockVector<double> solution_delta(m_dofs_per_block);
  solution_delta = 0.0;
 
  if (m_parameters.m_type_nonlinear_solver == "Newton")
  {
  bool newton_success = false;
  newton_success = solve_nonlinear_timestep_newton(solution_delta);
  AssertThrow(newton_success,
  ExcMessage("No convergence in Newton-Raphson nonlinear solver!"));
  /*

if Newton-Raphson failed, use LBFGS solver

  if (!newton_success)
  {
  solution_delta = 0.0;
  solve_nonlinear_timestep_LBFGS(solution_delta, LBFGS_update_refine);
  }
  */
  }
  else if (m_parameters.m_type_nonlinear_solver == "BFGS")
  solve_nonlinear_timestep_BFGS(solution_delta);
  else if (m_parameters.m_type_nonlinear_solver == "LBFGS")
  solve_nonlinear_timestep_LBFGS(solution_delta, LBFGS_update_refine);
  else
  AssertThrow(false, ExcMessage("Nonlinear solver type not implemented"));
 
  if (m_parameters.m_refinement_strategy == "adaptive-refine")
  {
 
  if (adp_refine_iteration == m_parameters.m_max_adaptive_refine_times)
  {
  m_solution += solution_delta;
  break;
  }
 
  mesh_is_same = local_refine_and_solution_transfer(solution_delta,
  LBFGS_update_refine);
 
  if (mesh_is_same)
  {
  m_solution += solution_delta;
  break;
  }
  }
  else if (m_parameters.m_refinement_strategy == "pre-refine")
  {
  m_solution += solution_delta;
  break;
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Selected mesh refinement strategy not implemented!"));
  }
  } // for (; adp_refine_iteration < m_parameters.m_max_adaptive_refine_times; ++adp_refine_iteration)
 

AssertThrow(adp_refine_iteration < m_parameters.m_max_adaptive_refine_times, ExcMessage("Number of local adaptive mesh refinement exceeds allowed maximum times!"));

  update_history_field_step();

output vtk files every 10 steps if there are too many time steps if (m_time.get_timestep() % 10 == 0)

  output_results();
 
  double energy_functional_current = calculate_energy_functional();
  m_logfile << "\t\tEnergy functional (J) = " << std::fixed << std::setprecision(10) << std::scientific
  << energy_functional_current << std::endl;
 
  std::pair<double, double> energy_pair = calculate_total_strain_energy_and_crack_energy_dissipation();
  m_logfile << "\t\tTotal strain energy (J) = " << std::fixed << std::setprecision(10) << std::scientific
  << energy_pair.first << std::endl;
  m_logfile << "\t\tCrack energy dissipation (J) = " << std::fixed << std::setprecision(10) << std::scientific
  << energy_pair.second << std::endl;
 
  std::pair<double, std::array<double, 3>> time_energy;
  time_energy.first = m_time.current();
  time_energy.second[0] = energy_pair.first;
  time_energy.second[1] = energy_pair.second;
  time_energy.second[2] = energy_pair.first + energy_pair.second;
  m_history_energy.push_back(time_energy);
 
  int face_ID = m_parameters.m_reaction_force_face_id;
  if (face_ID >= 0)
  calculate_reaction_force(face_ID);
 
  write_history_data();
 
  m_time.increment(time_table);
  } // while(m_time.current() < m_time.end() + m_time.get_delta_t()*1.0e-6)
  }
  } // namespace PhaseField
 
 
  int main(int argc, char* argv[])
  {
 
  using namespace dealii;
 
  if (argc != 2)
  AssertThrow(false,
  ExcMessage("The number of arguments provided to the program has to be 2!"));
 
  const unsigned int dim = std::stoi(argv[1]);
  if (dim == 2 )
  {
  PhaseField::PhaseFieldMonolithicSolve<2> FEQ1Full("parameters.prm");
  FEQ1Full.run();
  }
  else if (dim == 3)
  {
  PhaseField::PhaseFieldMonolithicSolve<3> SphereInclusion3D("parameters.prm");
  SphereInclusion3D.run();
  }
  else
  {
  AssertThrow(false,
  ExcMessage("Dimension has to be either 2 or 3"));
  }
 
  return 0;
  }
main
*  *  int main(int argc, char **argv)
Definition ALE_Finite_Strain_Plasticity.h:10156
Generated by doxygen 1.9.8