LifeV
HolzapfelMaterialNonLinear.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 _HOLZAPFELMATERIAL_H_
39 #define _HOLZAPFELMATERIAL_H_
40 
41 #pragma GCC diagnostic ignored "-Wunused-variable"
42 #pragma GCC diagnostic ignored "-Wunused-parameter"
43 
44 
45 #include <lifev/structure/solver/anisotropic/StructuralAnisotropicConstitutiveLaw.hpp>
46 
47 #define PI 3.14159265359
48 
49 namespace LifeV
50 {
51 template <typename MeshType>
52 class HolzapfelMaterialNonLinear : public StructuralAnisotropicConstitutiveLaw<MeshType>
53 {
54  //!@name Type definitions
55  //@{
56 
57 public:
58  typedef StructuralAnisotropicConstitutiveLaw<MeshType> super;
59 
60  typedef typename super::data_Type data_Type;
61 
62  typedef typename super::vector_Type vector_Type;
63  typedef typename super::matrix_Type matrix_Type;
64 
65  typedef typename super::matrixPtr_Type matrixPtr_Type;
66  typedef typename super::vectorPtr_Type vectorPtr_Type;
67  typedef typename super::dataPtr_Type dataPtr_Type;
69 
73 
77 
80 
83 
84  typedef typename super::FESpace_Type FESpace_Type;
87 
88  //Vector for vector parameters
91 
92  typedef MatrixSmall<3, 3> matrixSmall_Type;
93 
96 
99 
100  // Typedefs for expression definitions
101  typedef typename super::tensorF_Type tensorF_Type;
103  typedef typename super::tensorC_Type tensorC_Type;
104  typedef typename super::minusT_Type minusT_Type;
107 
108  // Anisotropic typedefs
111  typedef typename super::stretch_Type stretch_Type;
113  //@}
114 
115 
116 
117  //! @name Constructor & Destructor
118  //@{
119 
121 
122  virtual ~HolzapfelMaterialNonLinear();
123 
124  //@}
125 
126 
127 
128  //!@name Methods
129  //@{
130 
131  //! Setup the created object of the class StructuralIsotropicConstitutiveLaw
132  /*!
133  \param dFespace: the FiniteElement Space
134  \param monolithicMap: the MapEpetra
135  \param offset: the offset parameter used assembling the matrices
136  */
137  void setup ( const FESpacePtr_Type& dFESpace,
138  const ETFESpacePtr_Type& dETFESpace,
139  const std::shared_ptr<const MapEpetra>& monolithicMap,
140  const UInt offset,const dataPtr_Type& dataMaterial);
141 
142 
143  //! Compute the Stiffness matrix in StructuralSolver::buildSystem()
144  /*!
145  \param dataMaterial the class with Material properties data
146  */
147  void computeLinearStiff ( dataPtr_Type& /*dataMaterial*/,
148  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
149  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/);
150 
151 
152  //! Updates the Jacobian matrix in StructualSolver::updateJacobian
153  /*!
154  \param disp: solution at the k-th iteration of NonLinearRichardson Method
155  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
156  the material coefficients (e.g. Young modulus, Poisson ratio..)
157  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
158  */
159  void updateJacobianMatrix ( const vector_Type& disp,
160  const dataPtr_Type& dataMaterial,
161  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
162  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
163  const displayerPtr_Type& displayer );
164 
165 
166  //! Updates the nonlinear terms in the Jacobian matrix in StructualSolver::updateJacobian
167  /*!
168  \param stiff: stiffness matrix provided from outside
169  \param disp: solution at the k-th iteration of NonLinearRichardson Method
170  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
171  the material coefficients (e.g. Young modulus, Poisson ratio..)
172  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
173  */
175  const vector_Type& disp,
176  const dataPtr_Type& dataMaterial,
177  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
178  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/,
179  const displayerPtr_Type& displayer );
180 
181 
182  //! Interface method to compute the new Stiffness matrix in StructuralSolver::evalResidual and in
183  //! StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
184  /*!
185  \param sol: the solution vector
186  \param factor: scaling factor used in FSI
187  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the
188  material coefficients (e.g. Young modulus, Poisson ratio..)
189  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
190  */
191  void computeStiffness ( const vector_Type& disp,
192  const UInt iter,
193  Real factor,
194  const dataPtr_Type& dataMaterial,
195  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
196  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/,
197  const displayerPtr_Type& displayer );
198 
199 
201  const dataPtr_Type& dataMaterial,
202  const displayerPtr_Type& displayer ){};
203 
204  //! Computes the new Stiffness vector for Neo-Hookean and Holzapfel materials in StructuralSolver
205  //! given a certain displacement field.
206  //! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem
207  //! since the matrix is the expression of the matrix is the same.
208  /*!
209  \param sol: the solution vector
210  \param factor: scaling factor used in FSI
211  \param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
212  the material coefficients (e.g. Young modulus, Poisson ratio..)
213  \param displayer: a pointer to the Dysplaier member in the StructuralSolver class
214  */
215 
216  //! Computes the deformation gradient F, the cofactor matrix Cof(F),
217  //! the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
218  //! This function is used in StructuralIsotropicConstitutiveLaw::computeStiffness
219  /*!
220  \param dk_loc: the elemental displacement
221  */
222  void computeKinematicsVariables ( const VectorElemental& dk_loc )
223  {
224  }
225 
226  //! Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F), the trace of C Tr(C).
227  /*!
228  \param dk_loc: local displacement vector
229  */
230  //void computeStress( const vector_Type& sol );
231 
232  //! ShowMe method of the class (saved on a file the stiffness vector and the jacobian)
233  void showMe ( std::string const& fileNameVectStiff,
234  std::string const& fileNameJacobain );
235 
236  //! Compute the First Piola Kirchhoff Tensor
237  /*!
238  \param firstPiola Epetra_SerialDenseMatrix that has to be filled
239  \param tensorF Epetra_SerialDenseMatrix the deformation gradient
240  \param cofactorF Epetra_SerialDenseMatrix cofactor of F
241  \param invariants std::vector with the invariants of C and the detF
242  \param material UInt number to get the material parameteres form the VenantElasticData class
243  */
244  void computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
245  const Epetra_SerialDenseMatrix& tensorF,
246  const Epetra_SerialDenseMatrix& cofactorF,
247  const std::vector<Real>& invariants,
248  const UInt marker);
249 
250  //! Compute the First Piola Kirchhoff Tensor
251  /*!
252  \param disp the displacement field from which we compute the fisrt piola-Kirchhoff tensor
253  \param sigma_1 the first column of the Cauchy stress tensor
254  \param sigma_2 the second column of the Cauchy stress tensor
255  \param sigma_3 the third column of the Cauchy stress tensor
256  */
257  void computeCauchyStressTensor ( const vectorPtr_Type disp,
258  const QuadratureRule& evalQuad,
259  vectorPtr_Type sigma_1,
260  vectorPtr_Type sigma_2,
261  vectorPtr_Type sigma_3);
262 
263 
264  void setupFiberDirections( vectorFiberFunctionPtr_Type vectorOfFibers );
265 
266 
267  void evaluateActivationFibers( const vector_Type& displacement,
268  vector_Type& fourthInvariant){}
269  //@}
270 
271  //! @name Get Methods
272  //@{
273 
274  //! Get the Stiffness matrix
276  {
277  return super::M_jacobian;
278  }
279 
280 
281  //! Get the stiffness vector
283  {
284  return M_stiff;
285  }
286 
287  //! Specific for multimechanism
288  vectorPtr_Type const selectionCriterion( const UInt /*i*/ ) const
289  {
290  vectorPtr_Type empty( new vector_Type( this->M_dispFESpace->map() ) );
291  return empty;
292  }
293 
294 
295  vectorPtr_Type const activationDisplacement( const UInt /*i*/ ) const
296  {
297  vectorPtr_Type empty( new vector_Type( this->M_dispFESpace->map() ) );
298  return empty;
299  }
300 
301  vectorPtr_Type const activatedUnitFiber( const UInt /*i*/ ) const
302  {
303  vectorPtr_Type empty( new vector_Type( this->M_dispFESpace->map() ) );
304  return empty;
305  }
306 
307  vectorPtr_Type const activatedDeterminant( const UInt /*i*/ ) const
308  {
309  vectorPtr_Type empty( new vector_Type( this->M_dispFESpace->map() ) );
310  return empty;
311  }
312 
313 
314 
315  void apply ( const vector_Type& sol, vector_Type& res,
316  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
317  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
318  const displayerPtr_Type displayer) ;
319 
320  //@}
321 
322 
323 
324 protected:
325 
326  //! construct the vectors for the parameters
327  /*!
328  \param VOID
329  \return VOID
330  */
331  // void setupVectorsParameters ( void );
332 
333 
334  //! Vector: stiffness non-linear
336 
337  //Create the indentity for F
339 
340 };
341 
342 
343 
344 
345 
346 template <typename MeshType>
348  super ( ),
349  M_stiff ( ),
350  M_identity ( )
351 {
352 }
353 
354 
355 
356 
357 
358 template <typename MeshType>
360 {}
361 
362 
363 
364 
365 
366 template <typename MeshType>
367 void
368 HolzapfelMaterialNonLinear<MeshType>::setup ( const FESpacePtr_Type& dFESpace,
369  const ETFESpacePtr_Type& dETFESpace,
370  const std::shared_ptr<const MapEpetra>& monolithicMap,
371  const UInt offset,
372  const dataPtr_Type& dataMaterial)
373 {
374 
375  this->M_dataMaterial = dataMaterial;
376  this->M_dispFESpace = dFESpace;
377  this->M_dispETFESpace = dETFESpace;
378  this->M_localMap = monolithicMap;
379  this->M_offset = offset;
380 
381  M_stiff.reset (new vector_Type (*this->M_localMap) );
382 
383  M_identity (0, 0) = 1.0;
384  M_identity (0, 1) = 0.0;
385  M_identity (0, 2) = 0.0;
386  M_identity (1, 0) = 0.0;
387  M_identity (1, 1) = 1.0;
388  M_identity (1, 2) = 0.0;
389  M_identity (2, 0) = 0.0;
390  M_identity (2, 1) = 0.0;
391  M_identity (2, 2) = 1.0;
392 
393  this->M_epsilon = this->M_dataMaterial->smoothness();
394  //this->setupVectorsParameters( );
395 }
396 
397 
398 // template <typename MeshType>
399 // void
400 // HolzapfelMaterialNonLinear<MeshType>::setupVectorsParameters ( void )
401 // {
402 // // Paolo Tricerri: February, 20th
403 // // In each class, the name of the parameters has to inserted in the law
404 // // TODO: move the saving of the material parameters to more abstract objects
405 // // such that in the class of the material we do not need to call explicitly
406 // // the name of the parameter.
407 
408 // // Number of volume on the local part of the mesh
409 // UInt nbElements = this->M_dispFESpace->mesh()->numVolumes();
410 
411 // // Parameter alpha
412 // // 1. resize the vector in the first element of the vector.
413 // (* (this->M_vectorsParameters) ) [0].resize ( nbElements );
414 
415 // // Parameter gamma
416 // (* (this->M_vectorsParameters) ) [1].resize ( nbElements );
417 
418 // // Parameter bulk
419 // (* (this->M_vectorsParameters) ) [2].resize ( nbElements );
420 
421 // for (UInt i (0); i < nbElements; i++ )
422 // {
423 // // Extracting the marker
424 // UInt markerID = this->M_dispFESpace->mesh()->element ( i ).markerID();
425 
426 // Real alpha = this->M_dataMaterial->alpha ( markerID );
427 // Real gamma = this->M_dataMaterial->gamma ( markerID );
428 // Real bulk = this->M_dataMaterial->bulk ( markerID );
429 
430 // ( (* (this->M_vectorsParameters) ) [0]) [ i ] = alpha;
431 // ( (* (this->M_vectorsParameters) ) [1]) [ i ] = gamma;
432 // ( (* (this->M_vectorsParameters) ) [2]) [ i ] = bulk;
433 // }
434 // }
435 
436 
437 template <typename MeshType>
438 void
440 {
441  // Allocating the right space for the vector of fiber function
442  this->M_vectorOfFibers.reset( new vectorFiberFunction_Type( ) );
443 
444  // In this method, the vector of fiber functions has to be properly set in the main
445  // of the test. The functions are then projected on a FESpace for the integration
446 
447  // Number of volume on the local part of the mesh
448  UInt nbFamilies = (*vectorOfFibers).size();
449 
450  ASSERT( nbFamilies == this->M_dataMaterial->numberFibersFamilies(),
451  " The number of families set in the test is different from the one in the data" );
452 
453  this->M_vectorOfFibers->resize( nbFamilies );
454 
455  for( UInt k(0); k < nbFamilies; k++ )
456  {
457  ( *(this->M_vectorOfFibers) )[ k ] = ( *vectorOfFibers )[ k ];
458  }
459 
460  // Setting the vectors that will be used
461  this->M_vectorInterpolated.resize( nbFamilies );
462 
463  for( UInt k(0); k < nbFamilies; k++ )
464  {
465  this->M_vectorInterpolated[ k ].reset( new vector_Type(*this->M_localMap) );
466  this->M_dispFESpace->interpolate ( *( ( *(this->M_vectorOfFibers) )[ k ] ) ,
467  * ( ( this->M_vectorInterpolated )[ k ] ),
468  0.0 );
469  }
470 
471 }
472 
473 
474 template <typename MeshType>
475 void HolzapfelMaterialNonLinear<MeshType>::computeLinearStiff (dataPtr_Type& /*dataMaterial*/,
476  const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/,
477  const mapMarkerIndexesPtr_Type /*mapsMarkerIndexes*/)
478 {
479  //! Empty method for exponential material
480 }
481 
482 
483 
484 
485 
486 template <typename MeshType>
488  const dataPtr_Type& dataMaterial,
489  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
490  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
491  const displayerPtr_Type& displayer )
492 {
493 
494  this->M_jacobian.reset (new matrix_Type (*this->M_localMap) );
495 
496  displayer->leaderPrint (" \n*********************************\n ");
497  updateNonLinearJacobianTerms (this->M_jacobian, disp, this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
498  displayer->leaderPrint (" \n*********************************\n ");
499 
500 }
501 
502 template <typename MeshType>
504  const vector_Type& disp,
505  const dataPtr_Type& dataMaterial,
506  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
507  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
508  const displayerPtr_Type& displayer )
509 {
510  using namespace ExpressionAssembly;
511 
512  // In the following definitions, the critical template argument is MapEpetra
513  // When other maps will be available in LifeV, the Holzapfel class and its mother
514  // should have as template the MeshType and the MapType.
515 
516  // Definition of F
517  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
518 
519  // Definition of J
520  determinantF_Type J = ExpressionDefinitions::determinantF( F );
521 
522  //Definition of C
523  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
524 
525  // Definition of J^-(2/3) = det( C ) using the isochoric/volumetric decomposition
526  powerExpression_Type Jel = ExpressionDefinitions::powerExpression( J , (-2.0/3.0) );
527 
528  // Definition of F^-T
529  minusT_Type F_T = ExpressionDefinitions::minusT( F );
530 
531  // Update the heaviside function for the stretch of the fibers
532  * (jacobian) *= 0.0;
533 
534  displayer->leaderPrint (" Non-Linear S - updating non linear terms in the Jacobian Matrix (N fibers Holzapfel): \n");
535 
536  for( UInt i(0); i < this->M_vectorInterpolated.size(); i++ )
537  {
538 
539  displayer->leaderPrint (" ", i + 1,"-th fiber family \n" );
540  // Defining the expression for the i-th fiber
541  // Definitions of the quantities which depend on the fiber directions e.g. I_4^i
542  interpolatedValue_Type fiberIth = ExpressionDefinitions::interpolateFiber( this->M_dispETFESpace, *(this->M_vectorInterpolated[ i ] ) );
543 
544  // Definition of the tensor M = ithFiber \otimes ithFiber
545  // At the moment, it's automatic that the method constructs the expression M = ithFiber \otimes ithFiber
546  // For a more general case, the file ExpressionDefinitions.hpp should be changed
547  outerProduct_Type Mith = ExpressionDefinitions::fiberTensor( fiberIth );
548 
549  // Definition of the fourth invariant : I_4^i = C:Mith
550  stretch_Type IVith = ExpressionDefinitions::fiberStretch( C, Mith );
551 
552  // Definition of the fouth isochoric invariant : J^(-2.0/3.0) * I_4^i
553  isochoricStretch_Type IVithBar = ExpressionDefinitions::isochoricFourthInvariant( Jel, IVith );
554 
555  // first term:
556  // (-4.0/3.0) * aplha_i * J^(-2.0/3.0) * \bar{I_4} * ( \bar{I_4} - 1 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
557  // (epsilon / PI) * ( 1/ (1 + epsilon^2 * ( \bar{I_4} - 1 )^2 ) ) * ( F^-T : dF ) ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
558  integrate( elements ( this->M_dispETFESpace->mesh() ),
559  this->M_dispFESpace->qr(),
560  this->M_dispETFESpace,
561  this->M_dispETFESpace,
562  value( -4.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel *
563  IVithBar * ( IVithBar - value(1.0) ) *
564  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
565  derAtan( IVithBar - value(1.0), this->M_epsilon, ( 1.0 / PI ) ) *
566  dot( F_T, grad(phi_j) ) *
567  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
568  ) >> jacobian;
569 
570 
571  // second term
572  // 2.0 * aplha_i * J^(-4.0/3.0) * ( \bar{I_4} - 1 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
573  // (epsilon / PI) * ( 1/ (1 + epsilon^2 * ( \bar{I_4} - 1 )^2 ) ) * ( dF^T*F : M + F^T*dF:M ) ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
574  integrate( elements ( this->M_dispETFESpace->mesh() ),
575  this->M_dispFESpace->qr(),
576  this->M_dispETFESpace,
577  this->M_dispETFESpace,
578  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) *
579  ( IVithBar - value(1.0) ) * Jel * Jel *
580  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
581  derAtan( IVithBar - value(1.0), this->M_epsilon, ( 1.0 / PI ) ) *
582  dot( transpose( grad(phi_j) ) * F + transpose(F) * grad(phi_j) , Mith ) *
583  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
584  ) >> jacobian;
585 
586 
587  // third term
588  // ( -4.0/3.0 ) * alpha_i * J^(-2.0/3.0) * ( \bar{I_4} - 1 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
589  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( F^-T : dF) * ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
590  integrate( elements ( this->M_dispETFESpace->mesh() ),
591  this->M_dispFESpace->qr(),
592  this->M_dispETFESpace,
593  this->M_dispETFESpace,
594  value( -4.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) *
595  Jel * ( IVithBar - value(1.0) ) *
596  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
597  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
598  dot( F_T , grad(phi_j) ) *
599  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
600  ) >> jacobian;
601 
602  // fourth term
603  // ( -4.0/3.0 ) * alpha_i * J^(-2.0/3.0) * \bar{I_4} * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
604  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( F^-T : dF) * ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
605  integrate( elements ( this->M_dispETFESpace->mesh() ),
606  this->M_dispFESpace->qr(),
607  this->M_dispETFESpace,
608  this->M_dispETFESpace,
609  value( -4.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) *
610  Jel * IVithBar *
611  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
612  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
613  dot( F_T , grad(phi_j) ) *
614  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
615  ) >> jacobian;
616 
617 
618  // fifth term
619  // 2.0 * aplha_i * J^(-4.0/3.0) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
620  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( dF^T*F : M + F^T*dF:M ) ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
621  integrate( elements ( this->M_dispETFESpace->mesh() ),
622  this->M_dispFESpace->qr(),
623  this->M_dispETFESpace,
624  this->M_dispETFESpace,
625  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * Jel *
626  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
627  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
628  dot( transpose( grad(phi_j) ) * F + transpose(F) * grad(phi_j) , Mith ) *
629  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
630  ) >> jacobian;
631 
632  // sixth term
633  // (-8.0/3.0) * aplha_i * gamma_i * J^(-2.0/3.0) * \bar{I_4} * ( \bar{I_4} - 1.0 )^2 * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
634  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( F^-T:dF ) ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
635  integrate( elements ( this->M_dispETFESpace->mesh() ),
636  this->M_dispFESpace->qr(),
637  this->M_dispETFESpace,
638  this->M_dispETFESpace,
639  value( -8.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) *
640  Jel * IVithBar * ( IVithBar - value(1.0) ) * ( IVithBar - value(1.0) ) *
641  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
642  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
643  dot( F_T , grad(phi_j) ) *
644  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
645  ) >> jacobian;
646 
647 
648  // seventh term
649  // 4.0 * aplha_i * gamma_i * J^(-4.0/3.0) * ( \bar{I_4} - 1.0 )^2 * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
650  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( dF^T*F : M + F^T*dF:M ) ( ( F * M - (1.0/3.0) * I_4 * F^-T ) : d\phi )
651  integrate( elements ( this->M_dispETFESpace->mesh() ),
652  this->M_dispFESpace->qr(),
653  this->M_dispETFESpace,
654  this->M_dispETFESpace,
655  value( 4.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) *
656  Jel * Jel * ( IVithBar - value(1.0) ) * ( IVithBar - value(1.0) ) *
657  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
658  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
659  dot( transpose( grad(phi_j) ) * F + transpose(F) * grad(phi_j) , Mith ) *
660  dot( F * Mith - value(1.0/3.0) * IVith * F_T, grad(phi_i) )
661  ) >> jacobian;
662 
663  // tenth term
664  // 2.0 * aplha_i * J^(-2.0/3.0) * ( \bar{I_4} - 1.0 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
665  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( dF*M : d\phi )
666  integrate( elements ( this->M_dispETFESpace->mesh() ),
667  this->M_dispFESpace->qr(),
668  this->M_dispETFESpace,
669  this->M_dispETFESpace,
670  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
671  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
672  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
673  dot( grad(phi_j) * Mith, grad(phi_i) )
674  ) >> jacobian;
675 
676  // eleventh term
677  // (2.0/3.0) * aplha_i * \bar{I_4} * ( \bar{I_4} - 1.0 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
678  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( F^-T * dF^T * F^-T : d\phi )
679  integrate( elements ( this->M_dispETFESpace->mesh() ),
680  this->M_dispFESpace->qr(),
681  this->M_dispETFESpace,
682  this->M_dispETFESpace,
683  value( 2.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * IVithBar * ( IVithBar - value(1.0) ) *
684  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
685  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
686  dot( F_T * transpose( grad(phi_j) ) * F_T, grad(phi_i) )
687  ) >> jacobian;
688 
689 
690  // twelveth term
691  // (-2.0/3.0) * aplha_i * J^(-2.0/3.0) * ( \bar{I_4} - 1.0 ) * exp( gamma_i * ( \bar{I_4} - 1 )^2 ) *
692  // ( ( 1 / PI ) * atan(\epsilon(\bar{I_4} - 1)) + 1/2 ) * ( dF^T * F + F^T * dF : Mith ) ( F^-T : d \phi)
693  integrate( elements ( this->M_dispETFESpace->mesh() ),
694  this->M_dispFESpace->qr(),
695  this->M_dispETFESpace,
696  this->M_dispETFESpace,
697  value( -2.0/3.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
698  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
699  atan( IVithBar - value(1.0), this->M_epsilon, ( 1 / PI ), (1.0/2.0) ) *
700  dot( transpose( grad(phi_j) ) * F + transpose(F) * grad(phi_j) , Mith ) *
701  dot( F_T , grad( phi_i ) )
702  ) >> jacobian;
703 
704 
705 
706  } // closing loop on fibers
707 
708  jacobian->globalAssemble();
709 }
710 
711 
712 
713 
714 
715 template <typename MeshType>
717  const UInt /*iter*/,
718  Real /*factor*/,
719  const dataPtr_Type& dataMaterial,
720  const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
721  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
722  const displayerPtr_Type& displayer )
723 {
724 
725  using namespace ExpressionAssembly;
726 
727  M_stiff.reset (new vector_Type (*this->M_localMap) );
728  * (M_stiff) *= 0.0;
729 
730  displayer->leaderPrint (" \n*********************************\n ");
731  displayer->leaderPrint (" Non-Linear S- Computing the Holzapfel nonlinear stiffness vector (N fibers Holzapfel)");
732  displayer->leaderPrint (" \n*********************************\n ");
733 
734  // For anisotropic part of the Piola-Kirchhoff is assemble summing up the parts of the
735  // Piola-Kirchhoff using the fiber index
736 
737  // Here the fiber vector at the quadrature node is deduced using the method
738  // Interpolate value of the expression template.
739 
740  // Defining quantities which depend of the displacement field and not on the fiber direction
741 
742  // Definition of F
743  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, disp, this->M_offset, this->M_identity );
744 
745  // Definition of J
746  determinantF_Type J = ExpressionDefinitions::determinantF( F );
747 
748  //Definition of C
749  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
750 
751  // Definition of J^-(2/3) = det( C ) using the isochoric/volumetric decomposition
752  powerExpression_Type Jel = ExpressionDefinitions::powerExpression( J , (-2.0/3.0) );
753 
754  // Definition of F^-T
755  minusT_Type F_T = ExpressionDefinitions::minusT( F );
756 
757  for( UInt i(0); i < this->M_vectorInterpolated.size() ; i++ )
758  {
759  // As in other classes, the specialization of the MapType = MapEpetra makes this expression
760  // not always usable. When other maps will be available in LifeV, the class should be re-templated.
761 
762  // Defining the expression for the i-th fiber
763  // Definitions of the quantities which depend on the fiber directions e.g. I_4^i
764  interpolatedValue_Type fiberIth = ExpressionDefinitions::interpolateFiber( this->M_dispETFESpace, *(this->M_vectorInterpolated[ i ] ) );
765 
766  // Definition of the tensor M = ithFiber \otimes ithFiber
767  // At the moment, it's automatic that the method constructs the expression M = ithFiber \otimes ithFiber
768  // For a more general case, the file ExpressionDefinitions.hpp should be changed
769  outerProduct_Type Mith = ExpressionDefinitions::fiberTensor( fiberIth );
770 
771  // Definition of the fourth invariant : I_4^i = C:Mith
772  stretch_Type IVith = ExpressionDefinitions::fiberStretch( C, Mith );
773 
774  // Definition of the fouth isochoric invariant : J^(-2.0/3.0) * I_4^i
775  isochoricStretch_Type IVithBar = ExpressionDefinitions::isochoricFourthInvariant( Jel, IVith );
776 
777  // First term:
778  // 2 alpha_i J^(-2.0/3.0) ( \bar{I_4} - 1 ) exp( gamma_i * (\bar{I_4} - 1)^2 ) * F M : \grad phi_i
779  // where alpha_i and gamma_i are the fiber parameters and M is the 2nd order tensor of type f_i \otimes \ f_i
780  integrate ( elements ( this->M_dispETFESpace->mesh() ),
781  this->M_dispFESpace->qr(),
782  this->M_dispETFESpace,
783  atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
784  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
785  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
786  dot( F * Mith, grad( phi_i ) )
787  ) >> this->M_stiff;
788 
789  // Second term:
790  // 2 alpha_i J^(-2.0/3.0) ( \bar{I_4} - 1 ) exp( gamma_i * (\bar{I_4} - 1)^2 ) * ( 1.0/3.0 * I_4 ) F^-T : \grad phi_i
791  // where alpha_i and gamma_i are the fiber parameters and M is the 2nd order tensor of type f_i \otimes \ f_i
792  integrate ( elements ( this->M_dispETFESpace->mesh() ),
793  this->M_dispFESpace->qr(),
794  this->M_dispETFESpace,
795  atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
796  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
797  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
798  value( -1.0/3.0 ) * IVith *
799  dot( F_T , grad( phi_i ) )
800  ) >> this->M_stiff;
801 
802 
803  }
804 
805  this->M_stiff->globalAssemble();
806 }
807 
808 
809 template <typename MeshType>
810 void HolzapfelMaterialNonLinear<MeshType>::showMe ( std::string const& fileNameStiff,
811  std::string const& fileNameJacobian )
812 {
813  this->M_stiff->spy (fileNameStiff);
815 }
816 
817 
818 template <typename MeshType>
819 void HolzapfelMaterialNonLinear<MeshType>::apply ( const vector_Type& sol,
820  vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
821  const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
822  const displayerPtr_Type displayer)
823 {
824  computeStiffness (sol, 0, 0, this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
825  res += *M_stiff;
826 }
827 
828 
829 template <typename MeshType>
830 void HolzapfelMaterialNonLinear<MeshType>::computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
831  const Epetra_SerialDenseMatrix& tensorF,
832  const Epetra_SerialDenseMatrix& cofactorF,
833  const std::vector<Real>& invariants,
834  const UInt marker)
835 {
836 
837  // Still Need to Define P
838 
839 }
840 
841 template <typename MeshType>
843  const QuadratureRule& evalQuad,
844  vectorPtr_Type sigma_1,
845  vectorPtr_Type sigma_2,
846  vectorPtr_Type sigma_3)
847 {
848 
849  using namespace ExpressionAssembly;
850 
851  // Definition of F
852  tensorF_Type F = ExpressionDefinitions::deformationGradient( this->M_dispETFESpace, *disp, this->M_offset, this->M_identity );
853 
854  // Definition of J
855  determinantF_Type J = ExpressionDefinitions::determinantF( F );
856 
857  //Definition of C
858  tensorC_Type C = ExpressionDefinitions::tensorC( transpose(F), F );
859 
860  // Definition of J^-(2/3) = det( C ) using the isochoric/volumetric decomposition
861  powerExpression_Type Jel = ExpressionDefinitions::powerExpression( J , (-2.0/3.0) );
862 
863  // Definition of F^-T
864  minusT_Type F_T = ExpressionDefinitions::minusT( F );
865 
866  for( UInt i(0); i < this->M_vectorInterpolated.size() ; i++ )
867  {
868  // As in other classes, the specialization of the MapType = MapEpetra makes this expression
869  // not always usable. When other maps will be available in LifeV, the class should be re-templated.
870 
871  // Defining the expression for the i-th fiber
872  // Definitions of the quantities which depend on the fiber directions e.g. I_4^i
873  interpolatedValue_Type fiberIth = ExpressionDefinitions::interpolateFiber( this->M_dispETFESpace, *(this->M_vectorInterpolated[ i ] ) );
874 
875  // Definition of the tensor M = ithFiber \otimes ithFiber
876  // At the moment, it's automatic that the method constructs the expression M = ithFiber \otimes ithFiber
877  // For a more general case, the file ExpressionDefinitions.hpp should be changed
878  outerProduct_Type Mith = ExpressionDefinitions::fiberTensor( fiberIth );
879 
880  // Definition of the fourth invariant : I_4^i = C:Mith
881  stretch_Type IVith = ExpressionDefinitions::fiberStretch( C, Mith );
882 
883  // Definition of the fouth isochoric invariant : J^(-2.0/3.0) * I_4^i
884  isochoricStretch_Type IVithBar = ExpressionDefinitions::isochoricFourthInvariant( Jel, IVith );
885 
886  // First term:
887  // 2 alpha_i J^(-2.0/3.0) ( \bar{I_4} - 1 ) exp( gamma_i * (\bar{I_4} - 1)^2 ) * F M : \grad phi_i
888  // where alpha_i and gamma_i are the fiber parameters and M is the 2nd order tensor of type f_i \otimes \ f_i
889 
890  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
891  evalQuad,
892  this->M_dispETFESpace,
893  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
894  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
895  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
896  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
897  F * Mith ) *
898  transpose(F) , 0 ), phi_i )
899  ) >> sigma_1;
900 
901  // Second term:
902  // 2 alpha_i J^(-2.0/3.0) ( \bar{I_4} - 1 ) exp( gamma_i * (\bar{I_4} - 1)^2 ) * ( 1.0/3.0 * I_4 ) F^-T : \grad phi_i
903  // where alpha_i and gamma_i are the fiber parameters and M is the 2nd order tensor of type f_i \otimes \ f_i
904  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
905  evalQuad,
906  this->M_dispETFESpace,
907  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
908  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
909  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
910  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
911  value( -1.0/3.0 ) * IVith *
912  F_T ) *
913  transpose(F) , 0 ), phi_i )
914  ) >> sigma_1;
915  sigma_1->globalAssemble();
916 
917  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
918  evalQuad,
919  this->M_dispETFESpace,
920  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
921  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
922  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
923  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
924  F * Mith ) *
925  transpose(F) , 1 ), phi_i )
926  ) >> sigma_2;
927 
928  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
929  evalQuad,
930  this->M_dispETFESpace,
931  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
932  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
933  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
934  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
935  value( -1.0/3.0 ) * IVith *
936  F_T ) *
937  transpose(F) , 1 ), phi_i )
938  ) >> sigma_2;
939  sigma_2->globalAssemble();
940 
941  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
942  evalQuad,
943  this->M_dispETFESpace,
944  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
945  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
946  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
947  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
948  F * Mith ) *
949  transpose(F) , 2 ), phi_i )
950  ) >> sigma_3;
951 
952  evaluateNode( elements( this->M_dispETFESpace->mesh() ),
953  evalQuad,
954  this->M_dispETFESpace,
955  meas_K * dot( vectorFromMatrix( ( 1 / J ) *
956  ( atan( IVithBar - value(1.0) , this->M_epsilon, ( 1 / PI ), ( 1.0/2.0 ) ) *
957  value( 2.0 ) * value( this->M_dataMaterial->ithStiffnessFibers( i ) ) * Jel * ( IVithBar - value(1.0) ) *
958  exp( value( this->M_dataMaterial->ithNonlinearityFibers( i ) ) * ( IVithBar- value(1.0) ) * ( IVithBar- value(1.0) ) ) *
959  value( -1.0/3.0 ) * IVith *
960  F_T ) *
961  transpose(F) , 2 ), phi_i )
962  ) >> sigma_3;
963  sigma_3->globalAssemble();
964  }
965 
966 }
967 
968 
969 template <typename MeshType>
970 inline StructuralAnisotropicConstitutiveLaw<MeshType>* createHolzapfelMaterialNonLinear()
971 {
972  return new HolzapfelMaterialNonLinear<MeshType >();
973 }
974 
975 namespace
976 {
978 }
979 
980 } //Namespace LifeV
981 
982 #endif /* __EXPONENTIALMATERIAL_H */
void assignFunction(bcBase_Type &base)
Assign the function to the base of the BCHandler.
matrixPtr_Type const stiffMatrix() const
Get the Stiffness matrix.
void computeCauchyStressTensor(const vectorPtr_Type disp, const QuadratureRule &evalQuad, vectorPtr_Type sigma_1, vectorPtr_Type sigma_2, vectorPtr_Type sigma_3)
Compute the First Piola Kirchhoff Tensor.
void computeStiffness(const vector_Type &disp, const UInt iter, Real factor, const dataPtr_Type &dataMaterial, const mapMarkerVolumesPtr_Type, const mapMarkerIndexesPtr_Type, const displayerPtr_Type &displayer)
Interface method to compute the new Stiffness matrix in StructuralSolver::evalResidual and in Structu...
void evaluateActivationFibers(const vector_Type &displacement, vector_Type &fourthInvariant)
super::mapMarkerIndexesPtr_Type mapMarkerIndexesPtr_Type
super::mapMarkerVolumes_Type mapMarkerVolumes_Type
std::shared_ptr< vectorVolumes_Type > vectorVolumesPtr_Type
#define PI
super::interpolatedValue_Type interpolatedValue_Type
super::fiberFunctionPtr_Type fiberFunctionPtr_Type
super::vectorsParametersPtr_Type vectorsParametersPtr_Type
std::shared_ptr< vectorIndexes_Type > vectorIndexesPtr_Type
super::vectorsParameters_Type vectorsParameters_Type
mapMarkerIndexes_Type::const_iterator mapIteratorIndex_Type
std::vector< typename MeshType::element_Type * > vectorVolumes_Type
void computeReferenceConfigurations(const vector_Type &disp, const dataPtr_Type &dataMaterial, const displayerPtr_Type &displayer)
super::isochoricStretch_Type isochoricStretch_Type
super::powerExpression_Type powerExpression_Type
StructuralAnisotropicConstitutiveLaw< MeshType > super
void computeKinematicsVariables(const VectorElemental &dk_loc)
Computes the new Stiffness vector for Neo-Hookean and Holzapfel materials in StructuralSolver given a...
super::mapMarkerIndexes_Type mapMarkerIndexes_Type
void apply(const vector_Type &sol, vector_Type &res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const mapMarkerIndexesPtr_Type mapsMarkerIndexes, const displayerPtr_Type displayer)
#define ASSERT(X, A)
Definition: LifeAssert.hpp:90
void updateNonLinearJacobianTerms(matrixPtr_Type &jacobian, const vector_Type &disp, const dataPtr_Type &dataMaterial, const mapMarkerVolumesPtr_Type, const mapMarkerIndexesPtr_Type, const displayerPtr_Type &displayer)
Updates the nonlinear terms in the Jacobian matrix in StructualSolver::updateJacobian.
mapMarkerVolumes_Type::const_iterator mapIterator_Type
vectorPtr_Type const activatedUnitFiber(const UInt) const
Get the activation displacement.
const ExpressionPhiJ phi_j
Simple function to be used in the construction of an expression.
vectorPtr_Type const activatedDeterminant(const UInt) const
Get the activation displacement.
void showMe(std::string const &fileNameVectStiff, std::string const &fileNameJacobain)
Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F)...
void setup(const FESpacePtr_Type &dFESpace, const ETFESpacePtr_Type &dETFESpace, const std::shared_ptr< const MapEpetra > &monolithicMap, const UInt offset, const dataPtr_Type &dataMaterial)
Setup the created object of the class StructuralIsotropicConstitutiveLaw.
double Real
Generic real data.
Definition: LifeV.hpp:175
super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type
void computeLinearStiff(dataPtr_Type &, const mapMarkerVolumesPtr_Type, const mapMarkerIndexesPtr_Type)
Compute the Stiffness matrix in StructuralSolver::buildSystem()
super::vectorFiberFunction_Type vectorFiberFunction_Type
vectorPtr_Type const selectionCriterion(const UInt) const
Specific for multimechanism.
super::vectorFiberFunctionPtr_Type vectorFiberFunctionPtr_Type
void updateJacobianMatrix(const vector_Type &disp, const dataPtr_Type &dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const mapMarkerIndexesPtr_Type mapsMarkerIndexes, const displayerPtr_Type &displayer)
Updates the Jacobian matrix in StructualSolver::updateJacobian.
ExpressionScalar value(const Real &myValue)
Simple function to be used in the construction of an expression.
vectorPtr_Type M_stiff
construct the vectors for the parameters
const ExpressionMeas meas_K
Instance to be used in the expressions.
vectorPtr_Type const stiffVector() const
Get the stiffness vector.
void setupFiberDirections(vectorFiberFunctionPtr_Type vectorOfFibers)
vectorPtr_Type const activationDisplacement(const UInt) const
Get the activation displacement.
QuadratureRule - The basis class for storing and accessing quadrature rules.
StructuralAnisotropicConstitutiveLaw< MeshType > * createHolzapfelMaterialNonLinear()
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
void computeLocalFirstPiolaKirchhoffTensor(Epetra_SerialDenseMatrix &firstPiola, const Epetra_SerialDenseMatrix &tensorF, const Epetra_SerialDenseMatrix &cofactorF, const std::vector< Real > &invariants, const UInt marker)
Compute the First Piola Kirchhoff Tensor.
ExpressionDphiI grad(const ExpressionPhiI &)
Simple function to be used in the construction of an expression.