libMesh
nonlinear_neohooke_cc.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // \author Robert Weidlich
21 // \date Copyright 2012
22 
23 #include "nonlinear_neohooke_cc.h"
24 
25 void NonlinearNeoHookeCurrentConfig::init_for_qp(VectorValue<Gradient> & grad_u,
26  unsigned int qp)
27 {
28  this->current_qp = qp;
29  F.zero();
30  S.zero();
31 
32  {
33  RealTensor invF;
34  invF.zero();
35  for (unsigned int i = 0; i < 3; ++i)
36  for (unsigned int j = 0; j < 3; ++j)
37  invF(i, j) += libmesh_real(grad_u(i)(j));
38 
39  F.add(invF.inverse());
40 
41  libmesh_assert_greater (F.det(), -TOLERANCE);
42  }
43 
45  this->calculate_tangent();
46 
47  this->calculate_stress();
48 }
49 
50 
51 
53 {
54  Real mu = E / (2 * (1 + nu));
55  Real lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
56 
57  Real detF = F.det();
58 
59  C_mat.resize(6, 6);
60  for (unsigned int i = 0; i < 3; ++i) {
61  for (unsigned int j = 0; j < 3; ++j) {
62  if (i == j) {
63  C_mat(i, j) = 2 * mu + lambda;
64  C_mat(i + 3, j + 3) = mu - 0.5 * lambda * (detF * detF - 1);
65  } else {
66  C_mat(i, j) = lambda * detF * detF;
67  }
68  }
69  }
70 }
71 
72 
74 {
75  double mu = E / (2.0 * (1.0 + nu));
76  double lambda = E * nu / ((1 + nu) * (1 - 2 * nu));
77 
78  Real detF = F.det();
79  RealTensor Ft = F.transpose();
80 
81  RealTensor C = Ft * F;
82  RealTensor b = F * Ft;
83  RealTensor identity;
84  identity(0, 0) = 1.0; identity(1, 1) = 1.0; identity(2, 2) = 1.0;
85  RealTensor invC = C.inverse();
86 
87  S = 0.5 * lambda * (detF * detF - 1) * invC + mu * (identity - invC);
88 
89  tau = (F * S) * Ft;
90  sigma = 1/detF * tau;
91 }
92 
93 void NonlinearNeoHookeCurrentConfig::get_residual(DenseVector<Real> & residuum,
94  unsigned int & i)
95 {
96  B_L.resize(3, 6);
97  DenseVector<Real> sigma_voigt(6);
98 
99  this->build_b_0_mat(i, B_L);
100 
101  tensor_to_voigt(sigma, sigma_voigt);
102 
103  B_L.vector_mult(residuum, sigma_voigt);
104 }
105 
107  DenseVector<Real> & vec)
108 {
109  vec(0) = tensor(0, 0);
110  vec(1) = tensor(1, 1);
111  vec(2) = tensor(2, 2);
112  vec(3) = tensor(0, 1);
113  vec(4) = tensor(1, 2);
114  vec(5) = tensor(0, 2);
115 
116 }
117 
119  unsigned int & i,
120  unsigned int & j)
121 {
122  stiffness.resize(3, 3);
123 
124  double G_IK = (sigma * dphi[i][current_qp]) * dphi[j][current_qp];
125  stiffness(0, 0) += G_IK;
126  stiffness(1, 1) += G_IK;
127  stiffness(2, 2) += G_IK;
128 
129  B_L.resize(3, 6);
130  this->build_b_0_mat(i, B_L);
131  B_K.resize(3, 6);
132  this->build_b_0_mat(j, B_K);
133 
136  B_L *= 1/F.det();
137 
138  stiffness += B_L;
139 }
140 
141 void NonlinearNeoHookeCurrentConfig::build_b_0_mat(int i, DenseMatrix<Real> & b_0_mat)
142 {
143  for (unsigned int ii = 0; ii < 3; ++ii)
144  b_0_mat(ii, ii) = dphi[i][current_qp](ii);
145 
146  b_0_mat(0, 3) = dphi[i][current_qp](1);
147  b_0_mat(1, 3) = dphi[i][current_qp](0);
148  b_0_mat(1, 4) = dphi[i][current_qp](2);
149  b_0_mat(2, 4) = dphi[i][current_qp](1);
150  b_0_mat(0, 5) = dphi[i][current_qp](2);
151  b_0_mat(2, 5) = dphi[i][current_qp](0);
152 }
T libmesh_real(T a)
void init_for_qp(VectorValue< Gradient > &grad_u, unsigned int qp)
Initialize the class for the given displacement gradient at the specified quadrature point...
bool calculate_linearized_stiffness
Flag to indicate if it is necessary to calculate values for stiffness matrix during initialization...
void get_residual(DenseVector< Real > &residuum, unsigned int &i)
Return the residual vector for the current state.
void zero()
Set all entries of the tensor to 0.
Definition: type_tensor.h:1248
void build_b_0_mat(int i, DenseMatrix< Real > &b_l_mat)
RealTensorValue RealTensor
static const Real TOLERANCE
virtual void right_multiply(const DenseMatrixBase< T > &M2) libmesh_override
Performs the operation: (*this) <- (*this) * M3.
Definition: dense_matrix.C:210
void get_linearized_stiffness(DenseMatrix< Real > &stiffness, unsigned int &i, unsigned int &j)
Return the stiffness matrix for the current state.
static void tensor_to_voigt(const RealTensor &tensor, DenseVector< Real > &vec)
TypeTensor< T > inverse() const
Definition: type_tensor.h:1022
void right_multiply_transpose(const DenseMatrix< T > &A)
Right multiplies by the transpose of the matrix A.
Definition: dense_matrix.C:257
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void resize(const unsigned int new_m, const unsigned int new_n)
Resize the matrix.
Definition: dense_matrix.h:776
const std::vector< std::vector< RealGradient > > & dphi
void vector_mult(DenseVector< T > &dest, const DenseVector< T > &arg) const
Performs the matrix-vector multiplication, dest := (*this) * arg.
Definition: dense_matrix.C:382