LifeV
NeoHookean.cpp
Go to the documentation of this file.
1 #include <lifev/fsi_blocks/solver/NeoHookean.hpp>
2 
3 namespace LifeV
4 {
5 
6 NeoHookean::NeoHookean( const commPtr_Type& communicator ):
7  M_comm(communicator),
8  M_bulk(0.0),
9  M_mu(0.0),
10  M_offset ( 0 )
11 {
12  M_identity (0, 0) = 1.0;
13  M_identity (0, 1) = 0.0;
14  M_identity (0, 2) = 0.0;
15  M_identity (1, 0) = 0.0;
16  M_identity (1, 1) = 1.0;
17  M_identity (1, 2) = 0.0;
18  M_identity (2, 0) = 0.0;
19  M_identity (2, 1) = 0.0;
20  M_identity (2, 2) = 1.0;
21  if ( M_comm->MyPID() == 0 )
22  {
23  std::cout << "\nUsing NeoHookean model for the structure\n";
24  }
25 }
26 
27 NeoHookean::~NeoHookean()
28 {}
29 
30 void
31 NeoHookean::setCoefficients ( const Real density, const Real young, const Real poisson)
32 {
33  // Copying the parameters
34 
35  M_density = density;
36 
37  // Evaluation of the lame coefficients
38 
39  Real lambda = ( young * poisson ) / ( ( 1.0 + poisson ) * ( 1.0 - 2.0 * poisson ) );
40 
41  M_mu = young / ( 2.0 * ( 1.0 + poisson ) );
42 
43  M_bulk = ( 2.0 / 3.0 ) * M_mu + lambda;
44 
45 
46 }
47 
48 void
49 NeoHookean::setup( const meshPtr_Type& mesh, const std::string dOrder)
50 {
51  // Create FE spaces
52 
53  M_dOrder = dOrder;
54 
55  M_displacementFESpace.reset (new FESpace<mesh_Type, map_Type> (mesh, M_dOrder, 3, M_comm) );
56 
57  M_displacementFESpace_ETA.reset( new ETFESpace_displacement ( M_displacementFESpace->mesh(), &(M_displacementFESpace->refFE()), M_comm));
58 }
59 
60 void
61 NeoHookean::evaluate_residual( const vectorPtr_Type& solution, const Real& coefficient, const vectorPtr_Type& csi, vectorPtr_Type& residual )
62 {
63  // coefficient is: 1.0/ ( M_dt * M_dt * M_structureTimeAdvance->get_beta() )
64 
65  if ( M_comm->MyPID() == 0 )
66  {
67  std::cout << " S- Assembling residual NeoHookean structure\n";
68  }
69 
70  residual.reset (new vector_Type ( M_displacementFESpace->map() ) );
71  residual->zero();
72 
73  vectorPtr_Type solution_rep ( new vector_Type ( *solution, Repeated ) );
74 
75  vectorPtr_Type csi_rep ( new vector_Type ( *csi, Repeated ) );
76 
77  using namespace ExpressionAssembly;
78 
79  // Definition of F
80  tensorF_Type F = ExpressionDefinitions::deformationGradient( M_displacementFESpace_ETA, *solution_rep, M_offset, M_identity );
81 
82  // Definition of J
83  determinantF_Type J = ExpressionDefinitions::determinantF( F );
84 
85  // Definition of tensor C
86  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
87 
88  // Definition of F^-T
89  minusT_Type F_T = ExpressionDefinitions::minusT( F );
90 
91  // Definition of tr( C )
92  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
93 
94  integrate ( elements ( M_displacementFESpace_ETA->mesh() ),
95  M_displacementFESpace->qr(),
96  M_displacementFESpace_ETA,
97  value ( coefficient * M_density ) * dot( value( M_displacementFESpace_ETA, *solution_rep ), phi_i ) -
98  value ( M_density ) * dot ( value ( M_displacementFESpace_ETA, *csi_rep ), phi_i ) +
99  value ( M_mu ) * pow (J, -2.0 / 3.0) * (dot ( F - value (1.0 / 3.0) * I_C * F_T, grad (phi_i) ) )
100  + value (1.0 / 2.0) * value ( M_bulk ) * ( pow ( J , 2.0) - J + log (J) ) * dot ( F_T, grad (phi_i) )
101  ) >> residual;
102 
103  residual->globalAssemble();
104 }
105 
106 void
107 NeoHookean::update_jacobian(const vectorPtr_Type& solution, const Real& coefficient, matrixPtr_Type& jacobian )
108 {
109  using namespace ExpressionAssembly;
110 
111  jacobian->zero();
112 
113  vectorPtr_Type solution_rep ( new vector_Type ( *solution, Repeated ) );
114 
115  // Definition of F
116  tensorF_Type F = ExpressionDefinitions::deformationGradient( M_displacementFESpace_ETA, *solution_rep, M_offset, M_identity );
117 
118  // Definition of J
119  determinantF_Type J = ExpressionDefinitions::determinantF( F );
120 
121  // Definition of tensor C
122  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
123 
124  // Definition of F^-T
125  minusT_Type F_T = ExpressionDefinitions::minusT( F );
126 
127  // Definition of tr( C )
128  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
129 
130  integrate ( elements ( M_displacementFESpace_ETA->mesh() ) ,
131  M_displacementFESpace->qr(),
132  M_displacementFESpace_ETA,
133  M_displacementFESpace_ETA,
134  /* Inertia term */
135  value ( coefficient * M_density ) * dot( phi_i, phi_j )
136  /*Isochoric Part*/
137  /*Stiffness matrix : int { -2/3 * mu * J^(-2/3) *( F^-T : \nabla \delta ) ( F : \nabla \v ) } */
138  + value (-2.0 / 3.0) * value ( M_mu ) * pow (J, - (2.0 / 3.0) ) * dot ( F_T , grad (phi_j) ) * dot ( F , grad (phi_i) )
139  /*Stiffness matrix : int { 2/9 * mu * ( Ic_iso )( F^-T : \nabla \delta ) ( F^-T : \nabla \v ) } */
140  + value (2.0 / 9.0) * value ( M_mu ) * pow( J, (-2.0/3.0) ) * I_C * dot ( F_T , grad (phi_j) ) * dot ( F_T , grad (phi_i) )
141  /* Stiffness matrix : int { mu * J^(-2/3) (\nabla \delta : \nabla \v)} */
142  + value ( M_mu ) * pow (J, - (2.0 / 3.0) ) * dot ( grad (phi_j), grad (phi_i) )
143  /* Stiffness matrix : int { -2/3 * mu * J^(-2/3) ( F : \nabla \delta ) ( F^-T : \nabla \v ) } */
144  + value (-2.0 / 3.0) * value ( M_mu ) * pow (J, - (2.0 / 3.0) ) * dot ( F , grad (phi_j) ) * dot ( F_T , grad (phi_i) )
145  /* Stiffness matrix : int { 1/3 * mu * Ic_iso * (F^-T [\nabla \delta]^t F^-T ) : \nabla \v } */
146  + value (1.0 / 3.0) * value ( M_mu ) * pow( J, -2.0/3.0) * I_C * dot ( ( F_T * transpose (grad (phi_j) ) * F_T ), grad (phi_i) )
147  /*Volumetric Part*/
148  + value ( 1.0 / 2.0 ) * value ( M_bulk ) * ( value (2.0) *pow (J, 2.0) - J + value (1.0) ) * dot ( F_T, grad (phi_j) ) * dot ( F_T, grad (phi_i) )
149  + value ( - 1.0 / 2.0 ) * value ( M_bulk ) * ( pow (J, 2.0) - J + log (J) ) * dot ( F_T * transpose (grad (phi_j) ) * F_T, grad (phi_i) )
150  ) >> jacobian;
151 
152  jacobian->globalAssemble();
153 
154 }
155 
156 }