LifeV
SecondOrderExponentialMaterialNonLinear.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published
12  by the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  * @file
29  * @brief This file contains the definition for the St. Venant Kirchhoff linear material
30  *
31  * @version 1.0
32  * @date 29-07-2010
33  * @author Paolo Tricerri,
34  *
35  * @maintainer Paolo Tricerri <paolo.tricerri@epfl.ch>
36  */
37 
38 #ifndef _SECONDORDEREXPONENTIALMATERIAL_H_
39 #define _SECONDORDEREXPONENTIALMATERIAL_H_
40 
41 #pragma GCC diagnostic ignored "-Wunused-variable"
42 #pragma GCC diagnostic ignored "-Wunused-parameter"
43 
44 #include <lifev/structure/solver/isotropic/StructuralIsotropicConstitutiveLaw.hpp>
45 
46 namespace LifeV
47 {
48 template <typename MeshType>
49 class SecondOrderExponentialMaterialNonLinear : public StructuralIsotropicConstitutiveLaw<MeshType>
50 {
51  //!@name Type definitions
52  //@{
53 
54 public:
55  typedef StructuralIsotropicConstitutiveLaw<MeshType> super;
56 
57  typedef typename super::data_Type data_Type;
58 
59  typedef typename super::vector_Type vector_Type;
60  typedef typename super::matrix_Type matrix_Type;
61 
62  typedef typename super::matrixPtr_Type matrixPtr_Type;
63  typedef typename super::vectorPtr_Type vectorPtr_Type;
64  typedef typename super::dataPtr_Type dataPtr_Type;
65  typedef typename super::displayerPtr_Type displayerPtr_Type;
66 
67  typedef typename super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type;
68  typedef typename super::mapMarkerVolumes_Type mapMarkerVolumes_Type;
69  typedef typename mapMarkerVolumes_Type::const_iterator mapIterator_Type;
70 
71  typedef typename super::vectorVolumes_Type vectorVolumes_Type;
72  typedef std::shared_ptr<vectorVolumes_Type> vectorVolumesPtr_Type;
73 
74  typedef typename super::mapMarkerIndexesPtr_Type mapMarkerIndexesPtr_Type;
75  typedef typename super::mapMarkerIndexes_Type mapMarkerIndexes_Type;
76  typedef typename mapMarkerIndexes_Type::const_iterator mapIteratorIndex_Type;
77 
78  typedef std::vector<UInt> vectorIndexes_Type;
79  typedef std::shared_ptr<vectorIndexes_Type> vectorIndexesPtr_Type;
80 
81  typedef typename super::FESpacePtr_Type FESpacePtr_Type;
82  typedef typename super::ETFESpacePtr_Type ETFESpacePtr_Type;
83 
84  //Vector for vector parameters
85  typedef typename super::vectorsParameters_Type vectorsParameters_Type;
86  typedef typename super::vectorsParametersPtr_Type vectorsParametersPtr_Type;
87 
88  typedef MatrixSmall<3, 3> matrixSmall_Type;
89 
90  // Typedefs for expression definitions
91  typedef typename super::tensorF_Type tensorF_Type;
92  typedef typename super::determinantF_Type determinantF_Type;
93  typedef typename super::tensorC_Type tensorC_Type;
94  typedef typename super::minusT_Type minusT_Type;
95  typedef typename super::traceTensor_Type traceTensor_Type;
96  //@}
97 
98 
99 
100  //! @name Constructor & Destructor
101  //@{
102 
103  SecondOrderExponentialMaterialNonLinear();
104 
105  virtual ~SecondOrderExponentialMaterialNonLinear();
106 
107  //@}
108 
109 
110 
111  //!@name Methods
112  //@{
113 
114  //! Setup the created object of the class StructuralIsotropicConstitutiveLaw
115  /*!
116  \param dFespace: the FiniteElement Space
117  \param monolithicMap: the MapEpetra
118  \param offset: the offset parameter used assembling the matrices
119  */
120  void setup ( const FESpacePtr_Type& dFESpace,
121  const ETFESpacePtr_Type& dETFESpace,
122  const std::shared_ptr<const MapEpetra>& monolithicMap,
123  const UInt offset, const dataPtr_Type& dataMaterial);
124 
125 
126  //! Compute the Stiffness matrix in StructuralSolver::buildSystem()
127  /*!
128  \param dataMaterial the class with Material properties data
129  */
130  void computeLinearStiff ( dataPtr_Type& /*dataMaterial*/, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/, const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/ );
131 
132 
133  //! Updates the Jacobian matrix in StructualSolver::updateJacobian
134  /*!
135  \param disp: solution at the k-th iteration of NonLinearRichardson Method
136  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
137  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
138  */
139  void updateJacobianMatrix ( const vector_Type& disp,
140  const dataPtr_Type& dataMaterial,
141  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
142  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
143  const displayerPtr_Type& displayer );
144 
145 
146  //! Updates the nonlinear terms in the Jacobian matrix in StructualSolver::updateJacobian
147  /*!
148  \param stiff: stiffness matrix provided from outside
149  \param disp: solution at the k-th iteration of NonLinearRichardson Method
150  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
151  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
152  */
153  void updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
154  const vector_Type& disp,
155  const dataPtr_Type& dataMaterial,
156  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
157  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/,
158  const displayerPtr_Type& displayer );
159 
160 
161  //! Interface method to compute the new Stiffness matrix in StructuralSolver::evalResidual and in
162  //! StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
163  /*!
164  \param sol: the solution vector
165  \param factor: scaling factor used in FSI
166  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
167  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
168  */
169  void computeStiffness ( const vector_Type& disp, Real factor, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
170  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/, const displayerPtr_Type& displayer );
171 
172 
173  //! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in StructuralSolver given a certain displacement field.
174  //! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
175  /*!
176  \param sol: the solution vector
177  \param factor: scaling factor used in FSI
178  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
179  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
180  */
181  // void computeVector( const vector_Type& sol,
182  // Real factor,
183  // const dataPtr_Type& dataMaterial,
184  // const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
185  // const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
186  // const displayerPtr_Type& displayer );
187 
188 
189  //! Computes the deformation gradient F, the cofactor matrix Cof(F), the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
190  //! This function is used in StructuralIsotropicConstitutiveLaw::computeStiffness
191  /*!
192  \param dk_loc: the elemental displacement
193  */
194  void computeKinematicsVariables ( const VectorElemental& dk_loc );
195 
196  //! Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F), the trace of C Tr(C).
197  /*!
198  \param dk_loc: local displacement vector
199  */
200  //void computeStress( const vector_Type& sol );
201 
202  //! ShowMe method of the class (saved on a file the stiffness vector and the jacobian)
203  void showMe ( std::string const& fileNameVectStiff,
204  std::string const& fileNameJacobain );
205 
206  //! Compute the First Piola Kirchhoff Tensor
207  /*!
208  \param firstPiola Epetra_SerialDenseMatrix that has to be filled
209  \param tensorF Epetra_SerialDenseMatrix the deformation gradient
210  \param cofactorF Epetra_SerialDenseMatrix cofactor of F
211  \param invariants std::vector with the invariants of C and the detF
212  \param material UInt number to get the material parameteres form the VenantElasticData class
213  */
214  void computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
215  const Epetra_SerialDenseMatrix& tensorF,
216  const Epetra_SerialDenseMatrix& cofactorF,
217  const std::vector<Real>& invariants,
218  const UInt marker);
219 
220  //! Compute the First Piola Kirchhoff Tensor
221  /*!
222  \param disp the displacement field from which we compute the fisrt piola-Kirchhoff tensor
223  \param sigma_1 the first column of the Cauchy stress tensor
224  \param sigma_2 the second column of the Cauchy stress tensor
225  \param sigma_3 the third column of the Cauchy stress tensor
226  */
227  void computeCauchyStressTensor ( const vectorPtr_Type disp,
228  const QuadratureRule& evalQuad,
229  vectorPtr_Type sigma_1,
230  vectorPtr_Type sigma_2,
231  vectorPtr_Type sigma_3);
232 
233  //@}
234 
235  //! @name Get Methods
236  //@{
237 
238  //! Get the Stiffness matrix
239  matrixPtr_Type const stiffMatrix() const
240  {
241  return super::M_jacobian;
242  }
243 
244 
245  //! Get the stiffness vector
246  vectorPtr_Type const stiffVector() const
247  {
248  return M_stiff;
249  }
250 
251  void apply ( const vector_Type& sol, vector_Type& res,
252  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
253  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
254  const displayerPtr_Type displayer) ;
255 
256  //@}
257 
258 
259 
260 protected:
261 
262  //! construct the vectors for the parameters
263  /*!
264  \param VOID
265  \return VOID
266  */
267  void setupVectorsParameters ( void );
268 
269  //! Vector: stiffness non-linear
270  vectorPtr_Type M_stiff;
271 
272  //Create the indentity for F
273  matrixSmall_Type M_identity;
274 };
275 
276 
277 
278 
279 
280 template <typename MeshType>
281 SecondOrderExponentialMaterialNonLinear<MeshType>::SecondOrderExponentialMaterialNonLinear() :
282  super ( ),
283  M_stiff ( ),
284  M_identity ( )
285 {
286 }
287 
288 
289 
290 
291 
292 template <typename MeshType>
293 SecondOrderExponentialMaterialNonLinear<MeshType>::~SecondOrderExponentialMaterialNonLinear()
294 {}
295 
296 
297 
298 
299 
300 template <typename MeshType>
301 void
302 SecondOrderExponentialMaterialNonLinear<MeshType>::setup (const FESpacePtr_Type& dFESpace,
303  const ETFESpacePtr_Type& dETFESpace,
304  const std::shared_ptr<const MapEpetra>& monolithicMap,
305  const UInt offset,
306  const dataPtr_Type& dataMaterial)
307 {
308 
309  this->M_dataMaterial = dataMaterial;
310  // std::cout<<"I am setting up the Material"<<std::endl;
311 
312  this->M_dispFESpace = dFESpace;
313  this->M_dispETFESpace = dETFESpace;
314  this->M_localMap = monolithicMap;
315  this->M_offset = offset;
316 
317  M_stiff.reset ( new vector_Type (*this->M_localMap) );
318 
319  M_identity (0, 0) = 1.0;
320  M_identity (0, 1) = 0.0;
321  M_identity (0, 2) = 0.0;
322  M_identity (1, 0) = 0.0;
323  M_identity (1, 1) = 1.0;
324  M_identity (1, 2) = 0.0;
325  M_identity (2, 0) = 0.0;
326  M_identity (2, 1) = 0.0;
327  M_identity (2, 2) = 1.0;
328 
329  // The 3 is because the law uses three parameters (alpha, gamma, bulk).
330  // another way would be to set up the number of constitutive parameters of the law
331  // in the data file to get the right size. Note the comment below.
332  this->M_vectorsParameters.reset ( new vectorsParameters_Type ( 3 ) );
333 
334  this->setupVectorsParameters();
335 }
336 
337 
338 template <typename MeshType>
339 void
340 SecondOrderExponentialMaterialNonLinear<MeshType>::setupVectorsParameters ( void )
341 {
342  // Paolo Tricerri: February, 20th
343  // In each class, the name of the parameters has to inserted in the law
344  // TODO: move the saving of the material parameters to more abstract objects
345  // such that in the class of the material we do not need to call explicitly
346  // the name of the parameter.
347 
348  // Number of volume on the local part of the mesh
349  UInt nbElements = this->M_dispFESpace->mesh()->numVolumes();
350 
351  // Parameter alpha
352  // 1. resize the vector in the first element of the vector.
353  (* (this->M_vectorsParameters) ) [0].resize ( nbElements );
354 
355  // Parameter gamma
356  (* (this->M_vectorsParameters) ) [1].resize ( nbElements );
357 
358  // Parameter bulk
359  (* (this->M_vectorsParameters) ) [2].resize ( nbElements );
360 
361  for (UInt i (0); i < nbElements; i++ )
362  {
363  // Extracting the marker
364  UInt markerID = this->M_dispFESpace->mesh()->element ( i ).markerID();
365 
366  Real alpha = this->M_dataMaterial->alpha ( markerID );
367  Real gamma = this->M_dataMaterial->gamma ( markerID );
368  Real bulk = this->M_dataMaterial->bulk ( markerID );
369 
370  ( (* (this->M_vectorsParameters) ) [0]) [ i ] = alpha;
371  ( (* (this->M_vectorsParameters) ) [1]) [ i ] = gamma;
372  ( (* (this->M_vectorsParameters) ) [2]) [ i ] = bulk;
373  }
374 }
375 
376 
377 template <typename MeshType>
378 void SecondOrderExponentialMaterialNonLinear<MeshType>::computeLinearStiff (dataPtr_Type& /*dataMaterial*/, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/, const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/)
379 {
380  //! Empty method for second order exponential material
381 }
382 
383 template <typename MeshType>
384 void SecondOrderExponentialMaterialNonLinear<MeshType>::updateJacobianMatrix ( const vector_Type& disp,
385  const dataPtr_Type& dataMaterial,
386  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
387  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
388  const displayerPtr_Type& displayer )
389 {
390  this->M_jacobian.reset (new matrix_Type (*this->M_localMap) );
391 
392  displayer->leaderPrint (" \n************************************************\n ");
393  updateNonLinearJacobianTerms (this->M_jacobian, disp, this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
394  displayer->leaderPrint (" \n************************************************\n ");
395 }
396 
397 
398 
399 
400 
401 template <typename MeshType>
402 void SecondOrderExponentialMaterialNonLinear<MeshType>::updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
403  const vector_Type& disp,
404  const dataPtr_Type& dataMaterial,
405  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
406  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
407  const displayerPtr_Type& displayer )
408 {
409 
410  using namespace ExpressionAssembly;
411 
412  *jacobian *= 0.0;
413 
414  displayer->leaderPrint (" Non-Linear S- updating non linear terms in the Jacobian Matrix (Second Order Exponential)");
415 
416  //! Nonlinear part of jacobian
417  //! loop on volumes (i)
418 
419  // mapIterator_Type it;
420  // //mapIteratorIndex_Type itIndex;
421 
422  // vectorVolumesPtr_Type pointerListOfVolumes;
423  // vectorIndexesPtr_Type pointerListOfIndexes;
424 
425  // for( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++ )
426  // {
427 
428  // //Given the marker pointed by the iterator, let's extract the material parameters
429  // UInt marker = it->first;
430 
431  // // Debug
432  // // UInt markerIndex = itIndex->first;
433  // // ASSERT( marker == markerIndex, "The list of volumes is referring to a marker that is not the same as the marker of index!!!");
434 
435  // pointerListOfVolumes.reset( new vectorVolumes_Type(it->second) );
436  // pointerListOfIndexes.reset( new vectorIndexes_Type( (*mapsMarkerIndexes)[marker] ) );
437 
438  // Real bulk = dataMaterial->bulk(marker);
439  // Real alpha = dataMaterial->alpha(marker);
440  // Real gamma = dataMaterial->gamma(marker);
441 
442  // Definition of F
443  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
444 
445  // Definition of J
446  determinantF_Type J = ExpressionDefinitions::determinantF( F );
447 
448  // Definition of tensor C
449  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
450 
451  // Definition of F^-T
452  minusT_Type F_T = ExpressionDefinitions::minusT( F );
453 
454  // Definition of tr( C )
455  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
456 
457  //! VOLUMETRIC PART
458  //! 1. Stiffness matrix: int { 1/2 * bulk * ( 2 - 1/J + 1/J^2 ) * ( CofF : \nabla \delta ) (CofF : \nabla v) }
459  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
460  this->M_dispFESpace->qr(),
461  this->M_dispETFESpace,
462  this->M_dispETFESpace,
463  value ( 1.0 / 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( value (2.0) *pow (J, 2.0) - J + value (1.0) ) * dot ( F_T, grad (phi_j) ) * dot ( F_T, grad (phi_i) )
464  ) >> jacobian;
465 
466  //! 2. Stiffness matrix: int { 1/2 * bulk * ( 1/J- 1 - log(J)/J^2 ) * ( CofF1 [\nabla \delta]^t CofF ) : \nabla v }
467  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
468  this->M_dispFESpace->qr(),
469  this->M_dispETFESpace,
470  this->M_dispETFESpace,
471  value ( - 1.0 / 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( pow (J, 2.0) - J + log (J) ) * dot ( F_T * transpose (grad (phi_j) ) * F_T, grad (phi_i) )
472  ) >> jacobian;
473 
474 
475  //! ISOCHORIC PART
476  //! 1. Stiffness matrix : int { (- 4/3 * alpha * J^(-2/3) * exp( gamma*( Ic_iso - 3)^2 )*
477  //! * ( Ic_isoK + ( Ic_isoK - 3 ) * ( 2 * gamma * ( Ic_isoK - 3) +1 ) ) ) * ( F^-T : \nabla \delta ) ( F : \nabla \v ) }
478  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
479  this->M_dispFESpace->qr(),
480  this->M_dispETFESpace,
481  this->M_dispETFESpace,
482  value ( -4.0 / 3.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow ( J, -2.0 / 3.0 ) *
483  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 )* I_C - value (3.0) ) ) *
484  ( pow( J, -2.0/3.0 ) * I_C + ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) + value (1.0) ) ) *
485  ( dot ( F_T, grad (phi_j) ) * dot ( F, grad (phi_i) ) )
486  ) >> jacobian;
487 
488 
489  //! 2. Stiffness matrix : int { 4 * alpha * J^(-4/3) * exp( gamma*( Ic_iso - 3)^2 ) * ( 1.0 + 2 * gamma * (Ic_isoK - 3 )^2 ) *
490  //! ( F : \nabla \delta ) ( F : \nabla \v ) }
491  integrate ( elements ( this->M_dispETFESpace->mesh() ),
492  this->M_dispFESpace->qr(),
493  this->M_dispETFESpace,
494  this->M_dispETFESpace,
495  value ( 4.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow ( J, (-4.0 / 3.0) ) *
496  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
497  ( value (1.0) + value ( 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
498  ( dot ( F, grad (phi_j) ) * dot ( F, grad (phi_i) ) )
499  ) >> jacobian;
500 
501  //! 3. Stiffness matrix :
502  //!int { ( 4.0/9.0 * alpha * Ic_isoK * exp( gamma*( Ic_iso - 3)^2 ) * ( Ic_isoK + 2 * gamma * (Ic_isoK - 3)^2 * Ic_isoK + (Ic_isoK - 3) ) ) *
503  //! ( F^-T : \nabla \delta ) ( F^-T : \nabla \v )}
504  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
505  this->M_dispFESpace->qr(),
506  this->M_dispETFESpace,
507  this->M_dispETFESpace,
508  value (4.0 / 9.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow( J, -2.0/3.0 ) * I_C *
509  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) ) *
510  ( ( pow( J, -2.0/3.0) * I_C - value (3.0) ) * ( value (1.0) + value (2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * pow(J, -2.0/3.0) * I_C * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) + pow( J, -2.0/3.0 ) * I_C ) *
511  ( dot ( F_T, grad (phi_j) ) * dot ( F_T, grad (phi_i) ) )
512  ) >> jacobian;
513 
514 
515  //! 4. Stiffness matrix :
516  //! int { (-4.0/3.0 * alpha * J^(-2/3) * exp( gamma*( Ic_iso - 3)*( Ic_iso - 3) ) * ( Ic_isoK + 2*gamma*(Ic_isok - 3)Ic + 1) ) *
517  //! ( F : \nabla \delta ) ( F^-T : \nabla \v ) }
518  integrate ( elements ( this->M_dispETFESpace->mesh() ),
519  this->M_dispFESpace->qr(),
520  this->M_dispETFESpace,
521  this->M_dispETFESpace,
522  value (-4.0 / 3.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow ( J, (-2.0 / 3.0) ) *
523  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
524  ( pow( J, -2.0/3.0 ) * I_C + ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( value ( 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * pow( J, -2.0/3.0 ) * I_C + value (1.0) ) ) *
525  ( dot ( F, grad (phi_j) ) * dot ( F_T, grad (phi_i) ) )
526  ) >> jacobian;
527 
528  //! 5. Stiffness matrix : int { (2 * alpha * J^(-2/3) * exp( gamma*( Ic_iso - 3)*( Ic_iso - 3))*(Ic_iso - 3) ) * (\nabla \delta: \nabla \v)}
529  integrate ( elements ( this->M_dispETFESpace->mesh() ),
530  this->M_dispFESpace->qr(),
531  this->M_dispETFESpace,
532  this->M_dispETFESpace,
533  value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow ( J, (-2.0 / 3.0) ) * ( pow(J, -2.0/3.0) * I_C - value (3.0) ) *
534  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow(J, -2.0/3.0) * I_C - value (3.0) ) * ( pow(J, -2.0/3.0) * I_C - value (3.0) ) ) *
535  dot ( grad (phi_j), grad (phi_i) )
536  ) >> jacobian;
537 
538  //! 6. Stiffness matrix : int { 2.0/3.0 * alpha * Ic_iso * ( Ic_iso - 3) * exp(gamma*( Ic_iso - 3)*( Ic_iso - 3)) *
539  //! (CofF [\nabla \delta]^t CofF ) : \nabla \v }
540  integrate ( elements ( this->M_dispETFESpace->mesh() ),
541  this->M_dispFESpace->qr(),
542  this->M_dispETFESpace,
543  this->M_dispETFESpace,
544  value (2.0 / 3.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow( J, -2.0/3.0 ) * I_C * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) *
545  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
546  dot ( F_T * transpose (grad (phi_j) ) * F_T , grad (phi_i) )
547  ) >> jacobian;
548 
549 
550  // }
551 
552  jacobian->globalAssemble();
553 }
554 
555 
556 
557 
558 
559 template <typename MeshType>
560 void SecondOrderExponentialMaterialNonLinear<MeshType>::computeStiffness ( const vector_Type& disp,
561  Real /*factor*/,
562  const dataPtr_Type& dataMaterial,
563  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
564  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
565  const displayerPtr_Type& displayer )
566 {
567  using namespace ExpressionAssembly;
568 
569  this->M_stiff.reset (new vector_Type (*this->M_localMap) );
570  * (M_stiff) *= 0.0;
571 
572  displayer->leaderPrint (" \n*********************************\n ");
573  displayer->leaderPrint (" Non-Linear S- Computing the Second Order Exponential nonlinear stiffness vector ");
574  displayer->leaderPrint (" \n*********************************\n ");
575 
576 
577  // mapIterator_Type it;
578  // //mapIteratorIndex_Type itIndex;
579 
580  // vectorVolumesPtr_Type pointerListOfVolumes;
581  // vectorIndexesPtr_Type pointerListOfIndexes;
582 
583 
584  // for( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++)
585  // {
586 
587  // //Given the marker pointed by the iterator, let's extract the material parameters
588  // UInt marker = it->first;
589 
590  // // Debug
591  // // UInt markerIndex = itIndex->first;
592  // // ASSERT( marker == markerIndex, "The list of volumes is referring to a marker that is not the same as the marker of index!!!");
593 
594  // pointerListOfVolumes.reset( new vectorVolumes_Type(it->second) );
595  // pointerListOfIndexes.reset( new vectorIndexes_Type( (*mapsMarkerIndexes)[marker] ) );
596 
597 
598  // Real bulk = dataMaterial->bulk(marker);
599  // Real alpha = dataMaterial->alpha(marker);
600  // Real gamma = dataMaterial->gamma(marker);
601 
602  // Definition of F
603  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
604 
605  // Definition of J
606  determinantF_Type J = ExpressionDefinitions::determinantF( F );
607 
608  // Definition of tensor C
609  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
610 
611  // Definition of F^-T
612  minusT_Type F_T = ExpressionDefinitions::minusT( F );
613 
614  // Definition of tr( C )
615  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
616 
617  //! Stiffness for non-linear terms of the Neo-Hookean model
618  /*!
619  The results of the integrals are stored at each step into elvecK, until to build K matrix of the bilinear form
620  */
621  //! Volumetric part
622  /*!
623  Source term Pvol: int { bulk /2* (J1^2 - J1 + log(J1) ) * 1/J1 * (CofF1 : \nabla v) }
624  */
625  //! Stiffness for non-linear terms of the Neo-Hookean model
626  //! Volumetric part : int { bulk /2* (J1^2 - J1 + log(J1) ) * 1/J1 * (CofF1 : \nabla v) }
627  integrate ( elements ( this->M_dispETFESpace->mesh() ),
628  this->M_dispFESpace->qr(),
629  this->M_dispETFESpace,
630  value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * dot ( F_T, grad (phi_i) )
631  ) >> M_stiff;
632 
633  //! Isochoric part
634  /*!
635  Source term P1iso_Exp: int { 2 * alpha * ( Ic1_iso - 3 ) * exp(gamma *( Ic1_iso -3 )^2) *
636  ( J1^(-2/3)* (F1 : \nabla v) - 1/3 * (Ic1_iso / J1) * (CofF1 : \nabla v) ) }
637  */
638  integrate ( elements ( this->M_dispETFESpace->mesh() ),
639  this->M_dispFESpace->qr(),
640  this->M_dispETFESpace,
641  value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) *
642  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
643  pow ( J, (-2.0 / 3.0) ) * dot ( F - value (1.0 / 3.0) * I_C * F_T , grad (phi_i) )
644  ) >> M_stiff;
645 
646 
647 
648  // }
649 
650  this->M_stiff->globalAssemble();
651 }
652 
653 
654 
655 
656 
657 template <typename MeshType>
658 void SecondOrderExponentialMaterialNonLinear<MeshType>::computeKinematicsVariables ( const VectorElemental& dk_loc )
659 {
660 
661  // Real s;
662 
663  // //! loop on quadrature points (ig)
664  // for ( UInt ig = 0; ig < this->M_dispFESpace->fe().nbQuadPt(); ig++ )
665  // {
666  // //! loop on space coordinates (icoor)
667  // for ( UInt icoor = 0; icoor < nDimensions; icoor++ )
668  // {
669  // //! loop on space coordinates (jcoor)
670  // for ( UInt jcoor = 0; jcoor < nDimensions; jcoor++ )
671  // {
672  // s = 0.0;
673  // for ( UInt i = 0; i < this->M_dispFESpace->fe().nbFEDof(); i++ )
674  // {
675  // //! \grad u^k at a quadrature point
676  // s += this->M_dispFESpace->fe().phiDer( i, jcoor, ig ) * dk_loc[ i + icoor * this->M_dispFESpace->fe().nbFEDof() ];
677  // }
678  // //! gradient of displacement
679  // (*M_Fk)[ icoor ][ jcoor ][ig ] = s;
680  // }
681  // }
682  // }
683 
684  // //! loop on quadrature points (ig)
685  // for ( UInt ig = 0; ig < this->M_dispFESpace->fe().nbQuadPt(); ig++ )
686  // {
687  // //! loop on space coordinates (icoor)
688  // for ( UInt icoor = 0;icoor < nDimensions; icoor++ )
689  // {
690  // //! deformation gradient Fk
691  // (*M_Fk)[ icoor ][ icoor ][ ig ] += 1.0;
692  // }
693  // }
694 
695  // Real a,b,c,d,e,f,g,h,i;
696 
697  // for( UInt ig=0; ig< this->M_dispFESpace->fe().nbQuadPt(); ig++ )
698  // {
699  // a = (*M_Fk)[ 0 ][ 0 ][ ig ];
700  // b = (*M_Fk)[ 0 ][ 1 ][ ig ];
701  // c = (*M_Fk)[ 0 ][ 2 ][ ig ];
702  // d = (*M_Fk)[ 1 ][ 0 ][ ig ];
703  // e = (*M_Fk)[ 1 ][ 1 ][ ig ];
704  // f = (*M_Fk)[ 1 ][ 2 ][ ig ];
705  // g = (*M_Fk)[ 2 ][ 0 ][ ig ];
706  // h = (*M_Fk)[ 2 ][ 1 ][ ig ];
707  // i = (*M_Fk)[ 2 ][ 2 ][ ig ];
708 
709  // //! determinant of deformation gradient Fk
710  // (*M_Jack)[ig] = a*( e*i - f*h ) - b*( d*i - f*g ) + c*( d*h - e*g );
711 
712  // ASSERT_PRE((*M_Jack)[ig] > 0, "Negative Jacobian. Error!" );
713 
714  // (*M_CofFk)[ 0 ][ 0 ][ ig ] = ( e*i - f*h );
715  // (*M_CofFk)[ 0 ][ 1 ][ ig ] = - ( d*i - g*f );
716  // (*M_CofFk)[ 0 ][ 2 ][ ig ] = ( d*h - e*g );
717  // (*M_CofFk)[ 1 ][ 0 ][ ig ] = - ( b*i - c*h );
718  // (*M_CofFk)[ 1 ][ 1 ][ ig ] = ( a*i - c*g );
719  // (*M_CofFk)[ 1 ][ 2 ][ ig ] = - ( a*h - g*b );
720  // (*M_CofFk)[ 2 ][ 0 ][ ig ] = ( b*f - c*e );
721  // (*M_CofFk)[ 2 ][ 1 ][ ig ] = - ( a*f - c*d );
722  // (*M_CofFk)[ 2 ][ 2 ][ ig ] = ( a*e - d*b );
723  // }
724 
725  // //! loop on quadrature points
726  // for ( UInt ig = 0;ig < this->M_dispFESpace->fe().nbQuadPt(); ig++ )
727  // {
728  // s = 0.0;
729  // for ( UInt i = 0; i < nDimensions; i++)
730  // {
731  // for ( UInt j = 0; j < nDimensions; j++)
732  // {
733  // //! trace of C1 = (F1k^t F1k)
734  // s += (*M_Fk)[ i ][ j ][ ig ] * (*M_Fk)[ i ][ j ][ ig ];
735  // }
736  // }
737  // (*M_trCk)[ ig ] = s;
738  // }
739 
740  // for ( UInt ig = 0; ig < this->M_dispFESpace->fe().nbQuadPt(); ig++ )
741  // {
742  // //! trace of deviatoric C
743  // (*M_trCisok)[ ig ] = std::pow((*M_Jack)[ ig ], -2./3.) * (*M_trCk)[ ig ];
744  // }
745 }
746 
747 
748 
749 
750 
751 template <typename MeshType>
752 void SecondOrderExponentialMaterialNonLinear<MeshType>::showMe ( std::string const& fileNameStiff,
753  std::string const& fileNameJacobian )
754 {
755  this->M_stiff->spy (fileNameStiff);
757 }
758 
759 
760 template <typename MeshType>
761 void SecondOrderExponentialMaterialNonLinear<MeshType>::apply ( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const mapMarkerIndexesPtr_Type mapsMarkerIndexes,const displayerPtr_Type displayer )
762 {
763  computeStiffness (sol, 0, this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
764  res += *M_stiff;
765 }
766 
767 
768 template <typename MeshType>
769 void SecondOrderExponentialMaterialNonLinear<MeshType>::computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
770  const Epetra_SerialDenseMatrix& tensorF,
771  const Epetra_SerialDenseMatrix& cofactorF,
772  const std::vector<Real>& invariants,
773  const UInt marker)
774 {
775  //Get the material parameters
776  Real alpha = this->M_dataMaterial->alpha (marker);
777  Real gamma = this->M_dataMaterial->gamma (marker);
778  Real bulk = this->M_dataMaterial->bulk (marker);
779 
780 
781  //Computing the first term \alphaJ^{-2/3}[F-(1/3)tr(C)F^{-T}]exp(\gamma(tr(Ciso) - 3)
782  Epetra_SerialDenseMatrix firstTerm (tensorF);
783  Epetra_SerialDenseMatrix copyCofactorF (cofactorF);
784 
785  Real scale (0.0);
786  scale = -invariants[0] / 3.0;
787  copyCofactorF.Scale ( scale );
788  firstTerm += copyCofactorF;
789 
790  //Computation trace of the isochoric C
791  Real trCiso (0.0);
792  trCiso = std::pow (invariants[3], - (2.0 / 3.0) ) * invariants[0];
793 
794  Real coef ( 0.0 );
795  coef = 2.0 * alpha * std::pow (invariants[3], - (2.0 / 3.0) ) * ( trCiso - 3.0 ) * std::exp ( gamma * ( trCiso - 3 ) * ( trCiso - 3 ) );
796  firstTerm.Scale ( coef );
797 
798  //Computing the second term (volumetric part) J*(bulk/2)(J-1+(1/J)*ln(J))F^{-T}
799  Epetra_SerialDenseMatrix secondTerm (cofactorF);
800  Real secCoef (0);
801  secCoef = invariants[3] * (bulk / 2.0) * (invariants[3] - 1 + (1.0 / invariants[3]) * std::log (invariants[3]) );
802 
803  secondTerm.Scale ( secCoef );
804 
805  firstPiola += firstTerm;
806  firstPiola += secondTerm;
807 
808 }
809 
810 template <typename MeshType>
811 void SecondOrderExponentialMaterialNonLinear<MeshType>::computeCauchyStressTensor ( const vectorPtr_Type disp,
812  const QuadratureRule& evalQuad,
813  vectorPtr_Type sigma_1,
814  vectorPtr_Type sigma_2,
815  vectorPtr_Type sigma_3)
816 
817 {
818 
819  using namespace ExpressionAssembly;
820 
821  // Definition of F
822  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, *disp, this->M_offset, this->M_identity );
823 
824  // Definition of J
825  determinantF_Type J = ExpressionDefinitions::determinantF( F );
826 
827  // Definition of tensor C
828  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
829 
830  // Definition of F^-T
831  minusT_Type F_T = ExpressionDefinitions::minusT( F );
832 
833  // Definition of tr( C )
834  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
835 
836 
837  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
838  evalQuad,
839  this->M_dispETFESpace,
840  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
841  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + //vol
842  value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) *
843  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
844  pow ( J, (-2.0 / 3.0) ) * ( F - value (1.0 / 3.0) * I_C * F_T )
845  ) *
846  transpose( F ), 0 ), phi_i)
847  ) >> sigma_1;
848  sigma_1->globalAssemble();
849 
850  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
851  evalQuad,
852  this->M_dispETFESpace,
853  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
854  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + //vol
855  value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) *
856  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
857  pow ( J, (-2.0 / 3.0) ) * ( F - value (1.0 / 3.0) * I_C * F_T )
858  ) *
859  transpose( F ) , 1 ), phi_i)
860  ) >> sigma_2;
861  sigma_2->globalAssemble();
862 
863  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
864  evalQuad,
865  this->M_dispETFESpace,
866  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
867  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + //vol
868  value ( 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * ( pow( J, -2.0/3.0) * I_C - value (3.0) ) *
869  exp ( parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) * ( pow( J, -2.0/3.0 ) * I_C - value (3.0) ) ) *
870  pow ( J, (-2.0 / 3.0) ) * ( F - value (1.0 / 3.0) * I_C * F_T )
871  ) *
872  transpose( F ) , 2 ), phi_i)
873  ) >> sigma_3;
874  sigma_3->globalAssemble();
875 
876 }
877 
878 template <typename MeshType>
879 inline StructuralIsotropicConstitutiveLaw<MeshType>* createSecondOrderExponentialMaterialNonLinear()
880 {
881  return new SecondOrderExponentialMaterialNonLinear<MeshType >();
882 }
883 
884 namespace
885 {
887 }
888 
889 } //Namespace LifeV
890 
891 #endif /* __SECONDORDEREXPONENTIALMATERIAL_H */
void assignFunction(bcBase_Type &base)
Assign the function to the base of the BCHandler.
rightCauchyGreenTensor_Type tensorC(const ExpressionTranspose< deformationGradient_Type > tF, const deformationGradient_Type F)
determinantTensorF_Type determinantF(const deformationGradient_Type F)
const ExpressionPhiJ phi_j
Simple function to be used in the construction of an expression.
minusTransposedTensor_Type minusT(const deformationGradient_Type F)
traceTensor_Type traceTensor(const rightCauchyGreenTensor_Type C)
double Real
Generic real data.
Definition: LifeV.hpp:175
ExpressionScalar value(const Real &myValue)
Simple function to be used in the construction of an expression.
const ExpressionMeas meas_K
Instance to be used in the expressions.
QuadratureRule - The basis class for storing and accessing quadrature rules.
ExpressionDphiJ grad(const ExpressionPhiJ &)
Simple function to be used in the construction of an expression.
const ExpressionPhiI phi_i
Simple function to be used in the construction of an expression.
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191
ExpressionDphiI grad(const ExpressionPhiI &)
Simple function to be used in the construction of an expression.