40 #ifndef _EXPONENTIALMATERIAL_H_ 41 #define _EXPONENTIALMATERIAL_H_ 45 #include <lifev/structure/solver/isotropic/StructuralIsotropicConstitutiveLaw.hpp> 49 template <
typename MeshType>
50 class ExponentialMaterialNonLinear :
public StructuralIsotropicConstitutiveLaw<MeshType>
56 typedef StructuralIsotropicConstitutiveLaw<MeshType> super;
58 typedef typename super::data_Type data_Type;
60 typedef typename super::vector_Type vector_Type;
61 typedef typename super::matrix_Type matrix_Type;
63 typedef typename super::matrixPtr_Type matrixPtr_Type;
64 typedef typename super::vectorPtr_Type vectorPtr_Type;
65 typedef typename super::dataPtr_Type dataPtr_Type;
66 typedef typename super::displayerPtr_Type displayerPtr_Type;
68 typedef typename super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type;
69 typedef typename super::mapMarkerVolumes_Type mapMarkerVolumes_Type;
70 typedef typename mapMarkerVolumes_Type::const_iterator mapIterator_Type;
72 typedef typename super::mapMarkerIndexesPtr_Type mapMarkerIndexesPtr_Type;
73 typedef typename super::mapMarkerIndexes_Type mapMarkerIndexes_Type;
74 typedef typename mapMarkerIndexes_Type::const_iterator mapIteratorIndex_Type;
76 typedef std::vector<
typename MeshType::element_Type*> vectorVolumes_Type;
77 typedef std::shared_ptr<vectorVolumes_Type> vectorVolumesPtr_Type;
79 typedef std::vector<UInt> vectorIndexes_Type;
80 typedef std::shared_ptr<vectorIndexes_Type> vectorIndexesPtr_Type;
82 typedef typename super::FESpacePtr_Type FESpacePtr_Type;
83 typedef typename super::ETFESpacePtr_Type ETFESpacePtr_Type;
86 typedef typename super::vectorsParameters_Type vectorsParameters_Type;
87 typedef typename super::vectorsParametersPtr_Type vectorsParametersPtr_Type;
89 typedef MatrixSmall<3, 3> matrixSmall_Type;
92 typedef typename super::tensorF_Type tensorF_Type;
93 typedef typename super::determinantF_Type determinantF_Type;
94 typedef typename super::tensorC_Type tensorC_Type;
95 typedef typename super::minusT_Type minusT_Type;
96 typedef typename super::traceTensor_Type traceTensor_Type;
104 ExponentialMaterialNonLinear();
106 virtual ~ExponentialMaterialNonLinear();
121 void setup (
const FESpacePtr_Type& dFESpace,
122 const ETFESpacePtr_Type& dETFESpace,
123 const std::shared_ptr<
const MapEpetra>& monolithicMap,
124 const UInt offset,
const dataPtr_Type& dataMaterial);
130 void computeLinearStiff ( dataPtr_Type& ,
131 const mapMarkerVolumesPtr_Type ,
132 const mapMarkerIndexesPtr_Type );
142 void updateJacobianMatrix (
const vector_Type& disp,
143 const dataPtr_Type& dataMaterial,
144 const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
145 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
146 const displayerPtr_Type& displayer );
157 void updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
158 const vector_Type& disp,
159 const dataPtr_Type& dataMaterial,
160 const mapMarkerVolumesPtr_Type ,
161 const mapMarkerIndexesPtr_Type ,
162 const displayerPtr_Type& displayer );
174 void computeStiffness (
const vector_Type& disp,
176 const dataPtr_Type& dataMaterial,
177 const mapMarkerVolumesPtr_Type ,
178 const mapMarkerIndexesPtr_Type ,
179 const displayerPtr_Type& displayer );
205 void computeKinematicsVariables (
const VectorElemental& dk_loc )
216 void showMe ( std::string
const& fileNameVectStiff,
217 std::string
const& fileNameJacobain );
227 void computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
228 const Epetra_SerialDenseMatrix& tensorF,
229 const Epetra_SerialDenseMatrix& cofactorF,
230 const std::vector<Real>& invariants,
240 void computeCauchyStressTensor (
const vectorPtr_Type disp,
242 vectorPtr_Type sigma_1,
243 vectorPtr_Type sigma_2,
244 vectorPtr_Type sigma_3);
252 matrixPtr_Type
const stiffMatrix()
const 254 return super::M_jacobian;
259 vectorPtr_Type
const stiffVector()
const 264 void apply (
const vector_Type& sol, vector_Type& res,
265 const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
266 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
267 const displayerPtr_Type displayer) ;
280 void setupVectorsParameters (
void );
284 vectorPtr_Type M_stiff;
287 matrixSmall_Type M_identity;
295 template <
typename MeshType>
296 ExponentialMaterialNonLinear<MeshType>::ExponentialMaterialNonLinear() :
307 template <
typename MeshType>
308 ExponentialMaterialNonLinear<MeshType>::~ExponentialMaterialNonLinear()
315 template <
typename MeshType>
317 ExponentialMaterialNonLinear<MeshType>::setup (
const FESpacePtr_Type& dFESpace,
318 const ETFESpacePtr_Type& dETFESpace,
319 const std::shared_ptr<
const MapEpetra>& monolithicMap,
321 const dataPtr_Type& dataMaterial)
324 this->M_dataMaterial = dataMaterial;
325 this->M_dispFESpace = dFESpace;
326 this->M_dispETFESpace = dETFESpace;
327 this->M_localMap = monolithicMap;
328 this->M_offset = offset;
330 M_stiff.reset (
new vector_Type (*
this->M_localMap) );
332 M_identity (0, 0) = 1.0;
333 M_identity (0, 1) = 0.0;
334 M_identity (0, 2) = 0.0;
335 M_identity (1, 0) = 0.0;
336 M_identity (1, 1) = 1.0;
337 M_identity (1, 2) = 0.0;
338 M_identity (2, 0) = 0.0;
339 M_identity (2, 1) = 0.0;
340 M_identity (2, 2) = 1.0;
345 this->M_vectorsParameters.reset (
new vectorsParameters_Type ( 3 ) );
348 this->setupVectorsParameters( );
352 template <
typename MeshType>
354 ExponentialMaterialNonLinear<MeshType>::setupVectorsParameters (
void )
363 UInt nbElements =
this->M_dispFESpace->mesh()->numVolumes();
367 (* (
this->M_vectorsParameters) ) [0].resize ( nbElements );
370 (* (
this->M_vectorsParameters) ) [1].resize ( nbElements );
373 (* (
this->M_vectorsParameters) ) [2].resize ( nbElements );
375 for (
UInt i (0); i < nbElements; i++ )
378 UInt markerID =
this->M_dispFESpace->mesh()->element ( i ).markerID();
380 Real alpha =
this->M_dataMaterial->alpha ( markerID );
381 Real gamma =
this->M_dataMaterial->gamma ( markerID );
382 Real bulk =
this->M_dataMaterial->bulk ( markerID );
384 ( (* (
this->M_vectorsParameters) ) [0]) [ i ] = alpha;
385 ( (* (
this->M_vectorsParameters) ) [1]) [ i ] = gamma;
386 ( (* (
this->M_vectorsParameters) ) [2]) [ i ] = bulk;
391 template <
typename MeshType>
392 void ExponentialMaterialNonLinear<MeshType>::computeLinearStiff (dataPtr_Type& ,
393 const mapMarkerVolumesPtr_Type ,
394 const mapMarkerIndexesPtr_Type )
403 template <
typename MeshType>
404 void ExponentialMaterialNonLinear<MeshType>::updateJacobianMatrix (
const vector_Type& disp,
405 const dataPtr_Type& dataMaterial,
406 const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
407 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
408 const displayerPtr_Type& displayer )
411 this->M_jacobian.reset (
new matrix_Type (*
this->M_localMap) );
413 displayer->leaderPrint (
" \n*********************************\n ");
414 updateNonLinearJacobianTerms (
this->M_jacobian, disp,
this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
415 displayer->leaderPrint (
" \n*********************************\n ");
418 template <
typename MeshType>
419 void ExponentialMaterialNonLinear<MeshType>::updateNonLinearJacobianTerms ( matrixPtr_Type& jacobian,
420 const vector_Type& disp,
421 const dataPtr_Type& dataMaterial,
422 const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
423 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
424 const displayerPtr_Type& displayer )
428 displayer->leaderPrint (
" Non-Linear S- updating non linear terms in the Jacobian Matrix (Exponential)");
458 tensorF_Type F = ExpressionDefinitions::deformationGradient(
this->M_dispETFESpace, disp,
this->M_offset,
this->M_identity );
473 integrate ( elements (
this->M_dispETFESpace->mesh() ),
474 this->M_dispFESpace->qr(),
475 this->M_dispETFESpace,
476 this->M_dispETFESpace,
480 integrate ( elements (
this->M_dispETFESpace->mesh() ),
481 this->M_dispFESpace->qr(),
482 this->M_dispETFESpace,
483 this->M_dispETFESpace,
484 value ( -1.0 / 2.0 ) * parameter ( (* (
this->M_vectorsParameters) ) [2] ) * ( pow (J, 2.0) - J + log (J) ) * dot ( F_T * transpose (grad (phi_j) ) * F_T, grad (phi_i) )
490 integrate ( elements (
this->M_dispETFESpace->mesh() ),
491 this->M_dispFESpace->qr(),
492 this->M_dispETFESpace,
493 this->M_dispETFESpace,
494 value (-2.0 / 3.0
) * parameter ( (* (
this->M_vectorsParameters) ) [0]) *
495 pow (J, -2.0 / 3.0) * exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow( J, -(2.0/3.0) ) * I_C -
value (3.0
) ) ) *
496 (
value (1.0
) + parameter ( (* (
this->M_vectorsParameters) ) [1]) * pow( J, -(2.0/3.0) ) * I_C ) *
502 integrate ( elements (
this->M_dispETFESpace->mesh() ),
503 this->M_dispFESpace->qr(),
504 this->M_dispETFESpace,
505 this->M_dispETFESpace,
506 value (2.0
) * parameter ( (* (
this->M_vectorsParameters) ) [0] ) * parameter ( (* (
this->M_vectorsParameters) ) [1]) * pow (J, -4.0 / 3.0) * exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow( J, -(2.0/3.0) ) * I_C -
value (3.0
) ) ) *
512 integrate ( elements (
this->M_dispETFESpace->mesh() ) ,
513 this->M_dispFESpace->qr(),
514 this->M_dispETFESpace,
515 this->M_dispETFESpace,
516 value (2.0 / 9.0
) * parameter ( (* (
this->M_vectorsParameters) ) [0]) * ( pow( J, -(2.0/3.0) ) * I_C ) * exp ( parameter ( (* (
this->M_vectorsParameters) ) [1]) * ( pow( J, (-2.0/3.0) ) * I_C -
value (3.0
) ) ) * (
value (1.0
) + parameter ( (* (
this->M_vectorsParameters) ) [1]) * pow( J, (-2.0/3.0) ) * I_C ) *
522 integrate ( elements (
this->M_dispETFESpace->mesh() ),
523 this->M_dispFESpace->qr(),
524 this->M_dispETFESpace,
525 this->M_dispETFESpace,
526 value (-2.0 / 3.0
) * parameter ( (* (
this->M_vectorsParameters) ) [0]) * pow (J, -2.0 / 3.0) * exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow( J, (-2.0/3.0) ) * I_C -
value (3.0
) ) ) *
527 (
value (1.0
) + parameter ( (* (
this->M_vectorsParameters) ) [1]) * pow( J, (-2.0/3.0) ) * I_C ) *
533 integrate ( elements (
this->M_dispETFESpace->mesh() ),
534 this->M_dispFESpace->qr(),
535 this->M_dispETFESpace,
536 this->M_dispETFESpace,
537 parameter ( (* (
this->M_vectorsParameters) ) [0]) * pow (J, -2.0 / 3.0) *
538 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow( J, (-2.0/3.0) ) * I_C - value (3.0) ) ) *
539 dot ( grad (phi_j), grad (phi_i) )
545 integrate ( elements (
this->M_dispETFESpace->mesh() ),
546 this->M_dispFESpace->qr(),
547 this->M_dispETFESpace,
548 this->M_dispETFESpace,
549 value (1.0 / 3.0) * parameter ( (* (
this->M_vectorsParameters) ) [0]) * pow(J, (-2.0/3.0) ) * I_C *
550 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow(J, (-2.0/3.0) ) * I_C - value (3.0) ) ) *
551 dot ( F_T * transpose (grad (phi_j) ) * F_T , grad (phi_i) )
557 jacobian->globalAssemble();
564 template <
typename MeshType>
565 void ExponentialMaterialNonLinear<MeshType>::computeStiffness (
const vector_Type& disp,
567 const dataPtr_Type& dataMaterial,
568 const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
569 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
570 const displayerPtr_Type& displayer )
575 M_stiff.reset (
new vector_Type (*
this->M_localMap) );
578 displayer->leaderPrint (
" \n*********************************\n ");
579 displayer->leaderPrint (
" Non-Linear S- Computing the Exponential nonlinear stiffness vector ");
580 displayer->leaderPrint (
" \n*********************************\n ");
613 tensorF_Type F = ExpressionDefinitions::deformationGradient(
this->M_dispETFESpace, disp,
this->M_offset,
this->M_identity );
628 integrate ( elements (
this->M_dispETFESpace->mesh() ) ,
629 this->M_dispFESpace->qr(),
630 this->M_dispETFESpace,
631 value (1.0 / 2.0
) * parameter ( (* (
this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * dot ( F_T,
grad (phi_i) )
635 integrate ( elements (
this->M_dispETFESpace->mesh() ) ,
636 this->M_dispFESpace->qr(),
637 this->M_dispETFESpace,
638 parameter ( (* (
this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) *
639 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow (J, -2.0 / 3.0) * I_C -
value (3.0
) ) ) *
644 this->M_stiff->globalAssemble();
648 template <
typename MeshType>
657 template <
typename MeshType>
658 void ExponentialMaterialNonLinear<MeshType>::apply (
const vector_Type& sol,
659 vector_Type& res,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
660 const mapMarkerIndexesPtr_Type mapsMarkerIndexes,
const displayerPtr_Type displayer)
662 computeStiffness (sol, 0,
this->M_dataMaterial, mapsMarkerVolumes, mapsMarkerIndexes, displayer);
667 template <
typename MeshType>
668 void ExponentialMaterialNonLinear<MeshType>::computeLocalFirstPiolaKirchhoffTensor ( Epetra_SerialDenseMatrix& firstPiola,
669 const Epetra_SerialDenseMatrix& tensorF,
670 const Epetra_SerialDenseMatrix& cofactorF,
671 const std::vector<Real>& invariants,
676 Real alpha =
this->M_dataMaterial->alpha (marker);
677 Real gamma =
this->M_dataMaterial->gamma (marker);
678 Real bulk =
this->M_dataMaterial->bulk (marker);
682 Epetra_SerialDenseMatrix firstTerm (tensorF);
683 Epetra_SerialDenseMatrix copyCofactorF (cofactorF);
686 scale = -invariants[0] / 3.0;
687 copyCofactorF.Scale ( scale );
688 firstTerm += copyCofactorF;
692 trCiso = std::pow (invariants[3], - (2.0 / 3.0) ) * invariants[0];
695 coef = alpha * std::pow (invariants[3], - (2.0 / 3.0) ) * std::exp ( gamma * ( trCiso - 3 ) );
696 firstTerm.Scale ( coef );
699 Epetra_SerialDenseMatrix secondTerm (cofactorF);
701 secCoef = invariants[3] * (bulk / 2.0) * (invariants[3] - 1 + (1.0 / invariants[3]) * std::log (invariants[3]) );
703 secondTerm.Scale ( secCoef );
705 firstPiola += firstTerm;
706 firstPiola += secondTerm;
710 template <
typename MeshType>
711 void ExponentialMaterialNonLinear<MeshType>::computeCauchyStressTensor (
const vectorPtr_Type disp,
713 vectorPtr_Type sigma_1,
714 vectorPtr_Type sigma_2,
715 vectorPtr_Type sigma_3)
722 tensorF_Type F = ExpressionDefinitions::deformationGradient(
this->M_dispETFESpace, *disp,
this->M_offset,
this->M_identity );
736 evaluateNode( elements (
this->M_dispETFESpace->mesh() ),
738 this->M_dispETFESpace,
739 meas_K * dot ( vectorFromMatrix( ( 1 / J )*
740 (
value (1.0 / 2.0
) * parameter ( (* (
this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T +
741 parameter ( (* (
this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) *
742 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow (J, -2.0 / 3.0) * I_C -
value (3.0
) ) ) * ( F -
value(1.0/3.0
)* I_C* F_T ) ) *
743 transpose(F), 0 ),
phi_i)
746 sigma_1->globalAssemble();
747 evaluateNode( elements (
this->M_dispETFESpace->mesh() ),
749 this->M_dispETFESpace,
750 meas_K * dot ( vectorFromMatrix( ( 1 / J )*
751 (
value (1.0 / 2.0
) * parameter ( (* (
this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T +
752 parameter ( (* (
this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) *
753 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow (J, -2.0 / 3.0) * I_C -
value (3.0
) ) ) * ( F -
value(1.0/3.0
)* I_C* F_T ) ) *
754 transpose( F ) , 1 ),
phi_i)
757 sigma_2->globalAssemble();
759 evaluateNode( elements (
this->M_dispETFESpace->mesh() ),
761 this->M_dispETFESpace,
762 meas_K * dot ( vectorFromMatrix( ( 1 / J )*
763 (
value (1.0 / 2.0
) * parameter ( (* (
this->M_vectorsParameters) ) [2] ) * ( pow ( J , 2.0) - J + log (J) ) * F_T +
764 parameter ( (* (
this->M_vectorsParameters) ) [0] ) * pow (J, -2.0 / 3.0) *
765 exp ( parameter ( (* (
this->M_vectorsParameters) ) [1] ) * ( pow (J, -2.0 / 3.0) * I_C -
value (3.0
) ) ) * ( F -
value(1.0/3.0
)* I_C* F_T ) ) *
766 transpose( F ) , 2 ),
phi_i)
769 sigma_3->globalAssemble();
775 template <
typename MeshType>
776 inline StructuralIsotropicConstitutiveLaw<MeshType>* createExponentialMaterialNonLinear()
778 return new ExponentialMaterialNonLinear<MeshType >();
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)
static const LifeV::UInt elm_nodes_num[]
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.
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)
ExpressionDphiI grad(const ExpressionPhiI &)
Simple function to be used in the construction of an expression.