LifeV
NeoHookeanMaterialNonLinear.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  * @author Gianmarco Mengaldo
35  *
36  * @maintainer Paolo Tricerri <paolo.tricerri@epfl.ch>
37  * @contributor Gianmarco Mengaldo <gianmarco.mengaldo@gmail.com>
38  */
39 
40 #ifndef _NEOHOOKEANMATERIAL_H_
41 #define _NEOHOOKEANMATERIAL_H_
42 
43 
44 #include <lifev/structure/solver/isotropic/StructuralIsotropicConstitutiveLaw.hpp>
45 
46 namespace LifeV
47 {
48 
49 template <typename MeshType>
50 class NeoHookeanMaterialNonLinear :
51  public StructuralIsotropicConstitutiveLaw<MeshType>
52 {
53  //!@name Type definitions
54  //@{
55 
56 public:
57  typedef StructuralIsotropicConstitutiveLaw<MeshType> super;
58 
59  typedef StructuralConstitutiveLawData data_Type;
60 
61  typedef typename super::vector_Type vector_Type;
62  typedef typename super::matrix_Type matrix_Type;
63 
64  typedef typename super::matrixPtr_Type matrixPtr_Type;
65  typedef typename super::vectorPtr_Type vectorPtr_Type;
66  typedef typename super::dataPtr_Type dataPtr_Type;
67  typedef typename super::displayerPtr_Type displayerPtr_Type;
68 
69  typedef typename super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type;
70  typedef typename super::mapMarkerVolumes_Type mapMarkerVolumes_Type;
71  typedef typename mapMarkerVolumes_Type::const_iterator mapIterator_Type;
72 
73  typedef typename super::vectorVolumes_Type vectorVolumes_Type;
74  typedef std::shared_ptr<vectorVolumes_Type> vectorVolumesPtr_Type;
75 
76  typedef std::vector<UInt> vectorIndexes_Type;
77  typedef std::shared_ptr<vectorIndexes_Type> vectorIndexesPtr_Type;
78 
79  typedef typename super::mapMarkerIndexesPtr_Type mapMarkerIndexesPtr_Type;
80  typedef typename super::mapMarkerIndexes_Type mapMarkerIndexes_Type;
81  typedef typename mapMarkerIndexes_Type::const_iterator mapIteratorIndex_Type;
82 
83  typedef typename super::FESpacePtr_Type FESpacePtr_Type;
84  typedef typename super::ETFESpacePtr_Type ETFESpacePtr_Type;
85 
86  //Vector for vector parameters
87  typedef typename super::vectorsParameters_Type vectorsParameters_Type;
88  typedef typename super::vectorsParametersPtr_Type vectorsParametersPtr_Type;
89 
90  typedef MatrixSmall<3, 3> matrixSmall_Type;
91 
92  // Typedefs for expression definitions
93  typedef typename super::tensorF_Type tensorF_Type;
94  typedef typename super::determinantF_Type determinantF_Type;
95  typedef typename super::tensorC_Type tensorC_Type;
96  typedef typename super::minusT_Type minusT_Type;
97  typedef typename super::traceTensor_Type traceTensor_Type;
98  //@}
99 
100 
101 
102  //! @name Constructor & Destructor
103  //@{
104 
105  NeoHookeanMaterialNonLinear();
106 
107  virtual ~NeoHookeanMaterialNonLinear();
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*/,
131  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/ );
132 
133 
134  //! Updates the Jacobian matrix in StructualSolver::updateJacobian
135  /*!
136  \param disp: solution at the k-th iteration of NonLinearRichardson Method
137  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
138  the material coefficients (e.g. Young modulus, Poisson ratio..)
139  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
140  */
141  void updateJacobianMatrix ( const vector_Type& disp,
142  const dataPtr_Type& dataMaterial,
143  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
144  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
145  const displayerPtr_Type& displayer);
146 
147 
148  //! Updates the nonlinear terms in the Jacobian matrix in StructualSolver::updateJacobian
149  /*!
150  \param stiff: stiffness matrix provided from outside
151  \param disp: solution at the k-th iteration of NonLinearRichardson Method
152  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
153  the material coefficients (e.g. Young modulus, Poisson ratio..)
154  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
155  */
156  void updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
157  const vector_Type& disp,
158  const dataPtr_Type& dataMaterial,
159  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
160  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/,
161  const displayerPtr_Type& displayer);
162 
163  //! Interface method to compute the new Stiffness matrix in StructuralSolver::evalResidual and in
164  //! StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
165  /*!
166  \param sol: the solution vector
167  \param factor: scaling factor used in FSI
168  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
169  the material coefficients (e.g. Young modulus, Poisson ratio..)
170  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
171  */
172  void computeStiffness ( const vector_Type& disp, Real factor, const dataPtr_Type& dataMaterial,
173  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
174  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/,
175  const displayerPtr_Type& displayer );
176 
177  //! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in
178  //! StructuralSolver given a certain displacement field.
179  //! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem
180  //! since the matrix is the expression of the matrix is the same.
181  /*!
182  \param sol: the solution vector
183  \param factor: scaling factor used in FSI
184  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
185  the material coefficients (e.g. Young modulus, Poisson ratio..)
186  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
187  */
188  // void computeVector( const vector_Type& sol,
189  // Real factor,
190  // const dataPtr_Type& dataMaterial,
191  // const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
192  // const displayerPtr_Type& displayer );
193 
194 
195  //! Computes the deformation gradient F, the cofactor matrix Cof(F),
196  //! the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
197  //! This function is used in StructuralIsotropicConstitutiveLaw::computeStiffness
198  /*!
199  \param dk_loc: the elemental displacement
200  */
201  void computeKinematicsVariables ( const VectorElemental& dk_loc ) {}
202 
203  //! ShowMe method of the class (saved on a file the stiffness vector and the jacobian)
204  void showMe ( std::string const& fileNameVectStiff,
205  std::string const& fileNameJacobain);
206 
207 
208  //! Compute the First Piola Kirchhoff Tensor
209  /*!
210  \param firstPiola Epetra_SerialDenseMatrix that has to be filled
211  \param tensorF Epetra_SerialDenseMatrix the deformation gradient
212  \param cofactorF Epetra_SerialDenseMatrix cofactor of F
213  \param invariants std::vector with the invariants of C and the detF
214  \param material UInt number to get the material parameteres form the VenantElasticData class
215  */
216  void computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
217  const Epetra_SerialDenseMatrix& tensorF,
218  const Epetra_SerialDenseMatrix& cofactorF,
219  const std::vector<Real>& invariants,
220  const UInt marker);
221 
222  //! Compute the First Piola Kirchhoff Tensor
223  /*!
224  \param disp the displacement field from which we compute the fisrt piola-Kirchhoff tensor
225  \param sigma_1 the first column of the Cauchy stress tensor
226  \param sigma_2 the second column of the Cauchy stress tensor
227  \param sigma_3 the third column of the Cauchy stress tensor
228  */
229  void computeCauchyStressTensor ( const vectorPtr_Type disp,
230  const QuadratureRule& evalQuad,
231  vectorPtr_Type sigma_1,
232  vectorPtr_Type sigma_2,
233  vectorPtr_Type sigma_3);
234 
235  //@}
236 
237  //! @name Get Methods
238  //@{
239 
240  //! Get the Stiffness matrix
241  matrixPtr_Type const stiffMatrix() const
242  {
243  return super::M_jacobian;
244  }
245 
246  //! Get the stiffness vector
247  vectorPtr_Type const stiffVector() const
248  {
249  return M_stiff;
250  }
251 
252  void apply ( const vector_Type& sol, vector_Type& res,
253  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
254  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
255  const displayerPtr_Type displayer);
256 
257  //@}
258 
259 
260 
261 protected:
262 
263  //! construct the vectors for the parameters
264  /*!
265  \param VOID
266  \return VOID
267  */
268  void setupVectorsParameters ( void );
269 
270  //! Vector: stiffness non-linear
271  vectorPtr_Type M_stiff;
272 
273  //Create the indentity for F
274  matrixSmall_Type M_identity;
275 
276 };
277 
278 template <typename MeshType>
279 NeoHookeanMaterialNonLinear<MeshType>::NeoHookeanMaterialNonLinear() :
280  super ( ),
281  M_stiff ( ),
282  M_identity ( )
283 {
284 }
285 
286 
287 
288 
289 
290 template <typename MeshType>
291 NeoHookeanMaterialNonLinear<MeshType>::~NeoHookeanMaterialNonLinear()
292 {}
293 
294 
295 
296 
297 
298 template <typename MeshType>
299 void
300 NeoHookeanMaterialNonLinear<MeshType>::setup ( const FESpacePtr_Type& dFESpace,
301  const ETFESpacePtr_Type& dETFESpace,
302  const std::shared_ptr<const MapEpetra>& monolithicMap,
303  const UInt offset,
304  const dataPtr_Type& dataMaterial)
305 {
306  this->M_dataMaterial = dataMaterial;
307  this->M_dispFESpace = dFESpace;
308  this->M_dispETFESpace = dETFESpace;
309  this->M_localMap = monolithicMap;
310  this->M_offset = offset;
311  M_stiff.reset ( new vector_Type (*this->M_localMap) );
312 
313  M_identity (0, 0) = 1.0;
314  M_identity (0, 1) = 0.0;
315  M_identity (0, 2) = 0.0;
316  M_identity (1, 0) = 0.0;
317  M_identity (1, 1) = 1.0;
318  M_identity (1, 2) = 0.0;
319  M_identity (2, 0) = 0.0;
320  M_identity (2, 1) = 0.0;
321  M_identity (2, 2) = 1.0;
322 
323  // The 2 is because the law uses three parameters (mu, bulk).
324  // another way would be to set up the number of constitutive parameters of the law
325  // in the data file to get the right size. Note the comment below.
326  this->M_vectorsParameters.reset ( new vectorsParameters_Type ( 2 ) );
327 
328  this->setupVectorsParameters();
329 }
330 
331 template <typename MeshType>
332 void
333 NeoHookeanMaterialNonLinear<MeshType>::setupVectorsParameters ( void )
334 {
335  // Paolo Tricerri: February, 20th
336  // In each class, the name of the parameters has to inserted in the law
337  // TODO: move the saving of the material parameters to more abstract objects
338  // such that in the class of the material we do not need to call explicitly
339  // the name of the parameter.
340 
341  // Number of volume on the local part of the mesh
342  UInt nbElements = this->M_dispFESpace->mesh()->numVolumes();
343 
344  // Parameter mu
345  // 1. resize the vector in the first element of the vector.
346  (* (this->M_vectorsParameters) ) [0].resize ( nbElements );
347 
348  // Parameter bulk
349  (* (this->M_vectorsParameters) ) [1].resize ( nbElements );
350 
351  for (UInt i (0); i < nbElements; i++ )
352  {
353  // Extracting the marker
354  UInt markerID = this->M_dispFESpace->mesh()->element ( i ).markerID();
355 
356  Real mu = this->M_dataMaterial->mu ( markerID );
357  Real bulk = this->M_dataMaterial->bulk ( markerID );
358 
359  ( (* (this->M_vectorsParameters) ) [0]) [ i ] = mu;
360  ( (* (this->M_vectorsParameters) ) [1]) [ i ] = bulk;
361  }
362 }
363 
364 
365 template <typename MeshType>
366 void NeoHookeanMaterialNonLinear<MeshType>::computeLinearStiff (dataPtr_Type& /*dataMaterial*/,
367  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
368  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/)
369 {
370  //! Empty method for neo-hookean material
371 }
372 
373 
374 template <typename MeshType>
375 void NeoHookeanMaterialNonLinear<MeshType>::updateJacobianMatrix ( const vector_Type& disp,
376  const dataPtr_Type& dataMaterial,
377  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
378  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
379  const displayerPtr_Type& displayer )
380 {
381  this->M_jacobian.reset (new matrix_Type (*this->M_localMap) );
382 
383  displayer->leaderPrint (" \n*********************************\n ");
384  updateNonLinearJacobianTerms (this->M_jacobian, disp, this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
385  displayer->leaderPrint (" \n*********************************\n ");
386 }
387 
388 
389 template <typename MeshType>
390 void NeoHookeanMaterialNonLinear<MeshType>::updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
391  const vector_Type& disp,
392  const dataPtr_Type& dataMaterial,
393  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
394  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
395  const displayerPtr_Type& displayer )
396 {
397 
398  using namespace ExpressionAssembly;
399 
400  displayer->leaderPrint (" Non-Linear S- updating non linear terms in the Jacobian Matrix (Neo-Hookean)");
401 
402  * (jacobian) *= 0.0;
403 
404  // mapIterator_Type it;
405  // //mapIteratorIndex_Type itIndex;
406 
407  // vectorVolumesPtr_Type pointerListOfVolumes;
408  // vectorIndexesPtr_Type pointerListOfIndexes;
409 
410  // for( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++)
411  // {
412 
413  // //Given the marker pointed by the iterator, let's extract the material parameters
414  // UInt marker = it->first;
415 
416  // // Debug
417  // // UInt markerIndex = itIndex->first;
418  // // ASSERT( marker == markerIndex, "The list of volumes is referring to a marker that is not the same as the marker of index!!!");
419 
420  // pointerListOfVolumes.reset( new vectorVolumes_Type(it->second) );
421  // pointerListOfIndexes.reset( new vectorIndexes_Type( (*mapsMarkerIndexes)[marker]) );
422 
423  // Real mu = dataMaterial->mu(marker);
424  // Real bulk = dataMaterial->bulk(marker);
425 
426  // Definition of F
427  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
428 
429  // Definition of J
430  determinantF_Type J = ExpressionDefinitions::determinantF( F );
431 
432  // Definition of tensor C
433  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
434 
435  // Definition of F^-T
436  minusT_Type F_T = ExpressionDefinitions::minusT( F );
437 
438  // Definition of tr( C )
439  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
440 
441  //Assembling Volumetric Part
442  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
443  this->M_dispFESpace->qr(),
444  this->M_dispETFESpace,
445  this->M_dispETFESpace,
446  value ( 1.0 / 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( value (2.0) *pow (J, 2.0) - J + value (1.0) ) * dot ( F_T, grad (phi_j) ) * dot ( F_T, grad (phi_i) )
447  ) >> jacobian;
448 
449  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
450  this->M_dispFESpace->qr(),
451  this->M_dispETFESpace,
452  this->M_dispETFESpace,
453  value ( - 1.0 / 2.0 ) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow (J, 2.0) - J + log (J) ) * dot ( F_T * transpose (grad (phi_j) ) * F_T, grad (phi_i) )
454  ) >> jacobian;
455 
456 
457  //! ISOCHORIC PART
458  //! 1. Stiffness matrix : int { -2/3 * mu * J^(-2/3) *( F^-T : \nabla \delta ) ( F : \nabla \v ) }
459  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
460  this->M_dispFESpace->qr(),
461  this->M_dispETFESpace,
462  this->M_dispETFESpace,
463  value (-2.0 / 3.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, - (2.0 / 3.0) ) * dot ( F_T , grad (phi_j) ) * dot ( F , grad (phi_i) )
464  ) >> jacobian;
465 
466 
467  //! 2. Stiffness matrix : int { 2/9 * mu * ( Ic_iso )( F^-T : \nabla \delta ) ( F^-T : \nabla \v ) }
468  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
469  this->M_dispFESpace->qr(),
470  this->M_dispETFESpace,
471  this->M_dispETFESpace,
472  value (2.0 / 9.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) *
473  pow( J, (-2.0/3.0) ) * I_C * dot ( F_T , grad (phi_j) ) * dot ( F_T , grad (phi_i) )
474  ) >> jacobian;
475 
476  //! 3. Stiffness matrix : int { mu * J^(-2/3) (\nabla \delta : \nabla \v)}
477  integrate ( elements ( this->M_dispETFESpace->mesh() ) ,
478  this->M_dispFESpace->qr(),
479  this->M_dispETFESpace,
480  this->M_dispETFESpace,
481  parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, - (2.0 / 3.0) ) * dot ( grad (phi_j), grad (phi_i) )
482  ) >> jacobian;
483 
484  //! 4. Stiffness matrix : int { -2/3 * mu * J^(-2/3) ( F : \nabla \delta ) ( F^-T : \nabla \v ) }
485  integrate ( elements ( this->M_dispETFESpace->mesh() ),
486  this->M_dispFESpace->qr(),
487  this->M_dispETFESpace,
488  this->M_dispETFESpace,
489  value (-2.0 / 3.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) *
490  pow (J, - (2.0 / 3.0) ) * dot ( F , grad (phi_j) ) * dot ( F_T , grad (phi_i) )
491  ) >> jacobian;
492 
493 
494 
495  //! 5. Stiffness matrix : int { 1/3 * mu * Ic_iso * (F^-T [\nabla \delta]^t F^-T ) : \nabla \v }
496  integrate ( elements ( this->M_dispETFESpace->mesh() ),
497  this->M_dispFESpace->qr(),
498  this->M_dispETFESpace,
499  this->M_dispETFESpace,
500  value (1.0 / 3.0) * parameter ( (* (this->M_vectorsParameters) ) [0] ) *
501  pow( J, -2.0/3.0) * I_C * dot ( ( F_T * transpose (grad (phi_j) ) * F_T ), grad (phi_i) )
502  ) >> jacobian;
503 
504  // }
505 
506  jacobian->globalAssemble();
507 }
508 
509 
510 template <typename MeshType>
511 void NeoHookeanMaterialNonLinear<MeshType>::apply ( const vector_Type& sol, vector_Type& res,
512  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
513  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
514  const displayerPtr_Type displayer)
515 {
516  computeStiffness (sol, 0., this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
517  res += *M_stiff;
518 }
519 
520 
521 template <typename MeshType>
522 void NeoHookeanMaterialNonLinear<MeshType>::computeStiffness ( const vector_Type& disp,
523  Real /*factor*/,
524  const dataPtr_Type& dataMaterial,
525  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
526  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
527  const displayerPtr_Type& displayer )
528 {
529  using namespace ExpressionAssembly;
530 
531  this->M_stiff.reset (new vector_Type (*this->M_localMap) );
532 
533  displayer->leaderPrint (" \n******************************************************************\n ");
534  displayer->leaderPrint (" Non-Linear S- Computing the Neo-Hookean nonlinear stiffness vector" );
535  displayer->leaderPrint (" \n******************************************************************\n ");
536 
537  M_stiff.reset (new vector_Type (*this->M_localMap) );
538  * (M_stiff) *= 0.0;
539 
540  // mapIterator_Type it;
541  // //mapIteratorIndex_Type itIndex;
542 
543  // vectorVolumesPtr_Type pointerListOfVolumes;
544  // vectorIndexesPtr_Type pointerListOfIndexes;
545 
546  // for( it = (*mapsMarkerVolumes).begin(); it != (*mapsMarkerVolumes).end(); it++ )
547  // {
548 
549  // //Given the marker pointed by the iterator, let's extract the material parameters
550  // UInt marker = it->first;
551 
552  // // Debug
553  // // UInt markerIndex = itIndex->first;
554  // // ASSERT( marker == markerIndex, "The list of volumes is referring to a marker that is not the same as the marker of index!!!");
555 
556  // pointerListOfVolumes.reset( new vectorVolumes_Type(it->second) );
557  // pointerListOfIndexes.reset( new vectorIndexes_Type( (*mapsMarkerIndexes)[marker]) );
558 
559  // Real mu = dataMaterial->mu(marker);
560  // Real bulk = dataMaterial->bulk(marker);
561 
562  // Definition of F
563  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
564 
565  // Definition of J
566  determinantF_Type J = ExpressionDefinitions::determinantF( F );
567 
568  // Definition of tensor C
569  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
570 
571  // Definition of F^-T
572  minusT_Type F_T = ExpressionDefinitions::minusT( F );
573 
574  // Definition of tr( C )
575  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
576 
577  //Computation of the volumetric part
578  //
579  integrate ( elements ( this->M_dispETFESpace->mesh() ),
580  this->M_dispFESpace->qr(),
581  this->M_dispETFESpace,
582  value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow ( J , 2.0) - J + log (J) ) * dot ( F_T, grad (phi_i) )
583  ) >> M_stiff;
584 
585  //Computation of the isochoric part
586  integrate ( elements ( this->M_dispETFESpace->mesh() ),
587  this->M_dispFESpace->qr(),
588  this->M_dispETFESpace,
589  parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) * (dot ( F - value (1.0 / 3.0) * I_C * F_T, grad (phi_i) ) )
590  ) >> M_stiff;
591 
592  // }
593  this->M_stiff->globalAssemble();
594 }
595 
596 template <typename MeshType>
597 void NeoHookeanMaterialNonLinear<MeshType>::showMe ( std::string const& fileNameStiff,
598  std::string const& fileNameJacobian)
599 {
600  this->M_stiff->spy (fileNameStiff);
602 }
603 
604 template <typename MeshType>
605 void NeoHookeanMaterialNonLinear<MeshType>::computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
606  const Epetra_SerialDenseMatrix& tensorF,
607  const Epetra_SerialDenseMatrix& cofactorF,
608  const std::vector<Real>& invariants,
609  const UInt marker)
610 {
611 
612  //Get the material parameters
613  Real mu = this->M_dataMaterial->mu (marker);
614  Real bulk = this->M_dataMaterial->bulk (marker);
615 
616  //Computing the first term \muJ^{-2/3}[F-(1/3)tr(C)F^{-T}]
617  Epetra_SerialDenseMatrix firstTerm (tensorF);
618  Epetra_SerialDenseMatrix copyCofactorF (cofactorF);
619  Real scale ( 0.0 );
620  scale = -1 * (1.0 / 3.0) * invariants[0];
621  copyCofactorF.Scale ( scale );
622  firstTerm += copyCofactorF;
623 
624  Real coef ( 0.0 );
625  coef = mu * std::pow (invariants[3], -2.0 / 3.0);
626  firstTerm.Scale ( coef );
627 
628  //Computing the second term (volumetric part) J*(bulk/2)(J-1+(1/J)*ln(J))F^{-T}
629  Epetra_SerialDenseMatrix secondTerm (cofactorF);
630  Real sCoef (0);
631  sCoef = invariants[3] * (bulk / 2.0) * (invariants[3] - 1 + (1 / invariants[3]) * std::log (invariants[3]) );
632  secondTerm.Scale ( sCoef );
633 
634  firstPiola += firstTerm;
635  firstPiola += secondTerm;
636 }
637 
638 
639 template <typename MeshType>
640 void NeoHookeanMaterialNonLinear<MeshType>::computeCauchyStressTensor ( const vectorPtr_Type disp,
641  const QuadratureRule& evalQuad,
642  vectorPtr_Type sigma_1,
643  vectorPtr_Type sigma_2,
644  vectorPtr_Type sigma_3)
645 
646 {
647 
648  using namespace ExpressionAssembly;
649 
650  // Definition that are needed
651  // Definition of F
652  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, *disp, this->M_offset, this->M_identity );
653 
654  // Definition of J
655  determinantF_Type J = ExpressionDefinitions::determinantF( F );
656 
657  //Definition of C
658  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
659 
660  // Definition of F^-T
661  minusT_Type F_T = ExpressionDefinitions::minusT( F );
662 
663  // Definition I_C
664  traceTensor_Type I_C = ExpressionDefinitions::traceTensor( C );
665 
666  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
667  evalQuad,
668  this->M_dispETFESpace,
669  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
670  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + // volumetric
671  parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) * ( F - value (1.0 / 3.0) * I_C * F_T )
672  ) *
673  transpose(F), 0 ), phi_i)
674  ) >> sigma_1;
675  sigma_1->globalAssemble();
676 
677 
678  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
679  evalQuad,
680  this->M_dispETFESpace,
681  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
682  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + // volumetric
683  parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) * ( F - value (1.0 / 3.0) * I_C * F_T ) ) *
684  transpose( F ) , 1 ), phi_i)
685  ) >> sigma_2;
686  sigma_2->globalAssemble();
687 
688  evaluateNode( elements ( this->M_dispETFESpace->mesh() ),
689  evalQuad,
690  this->M_dispETFESpace,
691  meas_K * dot ( vectorFromMatrix( ( 1 / J )*
692  ( value (1.0 / 2.0) * parameter ( (* (this->M_vectorsParameters) ) [1] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T + // volumetric
693  parameter ( (* (this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) * ( F - value (1.0 / 3.0) * I_C * F_T )
694  ) *
695  transpose( F ) , 2 ), phi_i)
696  ) >> sigma_3;
697  sigma_3->globalAssemble();
698 
699 }
700 
701 
702 template <typename MeshType>
703 inline StructuralIsotropicConstitutiveLaw<MeshType>* createNeoHookeanMaterialNonLinear()
704 {
705  return new NeoHookeanMaterialNonLinear<MeshType >();
706 }
707 namespace
708 {
710 }
711 
712 } //Namespace LifeV
713 
714 #endif /* __NEOHOOKENANMATERIAL_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.
DataElasticStructure - Data container for solid problems with elastic structure.
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.