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