38 #ifndef _DISTRIBUTEDHOLZAPFELMATERIAL_H_ 39 #define _DISTRIBUTEDHOLZAPFELMATERIAL_H_ 41 #pragma GCC diagnostic ignored "-Wunused-variable" 42 #pragma GCC diagnostic ignored "-Wunused-parameter" 45 #include <lifev/structure/solver/anisotropic/StructuralAnisotropicConstitutiveLaw.hpp> 47 #define PI 3.14159265359
51 template <
typename MeshType>
58 typedef StructuralAnisotropicConstitutiveLaw<MeshType>
super;
146 const std::shared_ptr<
const MapEpetra>& monolithicMap,
251 const Epetra_SerialDenseMatrix& tensorF,
252 const Epetra_SerialDenseMatrix& cofactorF,
253 const std::vector<Real>& invariants,
283 return super::M_jacobian;
352 template <
typename MeshType>
364 template <
typename MeshType>
372 template <
typename MeshType>
376 const std::shared_ptr<
const MapEpetra>& monolithicMap,
381 this->M_dataMaterial = dataMaterial;
382 this->M_dispFESpace = dFESpace;
383 this->M_dispETFESpace = dETFESpace;
384 this->M_localMap = monolithicMap;
385 this->M_offset = offset;
389 M_identity (0, 0) = 1.0;
390 M_identity (0, 1) = 0.0;
391 M_identity (0, 2) = 0.0;
392 M_identity (1, 0) = 0.0;
393 M_identity (1, 1) = 1.0;
394 M_identity (1, 2) = 0.0;
395 M_identity (2, 0) = 0.0;
396 M_identity (2, 1) = 0.0;
397 M_identity (2, 2) = 1.0;
399 this->M_epsilon =
this->M_dataMaterial->smoothness();
443 template <
typename MeshType>
454 UInt nbFamilies = (*vectorOfFibers).size();
456 ASSERT( nbFamilies ==
this->M_dataMaterial->numberFibersFamilies(),
457 " The number of families set in the test is different from the one in the data" );
459 this->M_vectorOfFibers->resize( nbFamilies );
461 for(
UInt k(0); k < nbFamilies; k++ )
463 ( *(
this->M_vectorOfFibers) )[ k ] = ( *vectorOfFibers )[ k ];
467 this->M_vectorInterpolated.resize( nbFamilies );
469 for(
UInt k(0); k < nbFamilies; k++ )
471 this->M_vectorInterpolated[ k ].reset(
new vector_Type(*
this->M_localMap) );
472 this->M_dispFESpace->interpolate ( *( ( *(
this->M_vectorOfFibers) )[ k ] ) ,
473 * ( (
this->M_vectorInterpolated )[ k ] ),
480 template <
typename MeshType>
492 template <
typename MeshType>
500 this->M_jacobian.reset (
new matrix_Type (*
this->M_localMap) );
502 displayer->leaderPrint (
" \n*********************************\n ");
504 displayer->leaderPrint (
" \n*********************************\n ");
508 template <
typename MeshType>
523 tensorF_Type F = ExpressionDefinitions::deformationGradient(
this->M_dispETFESpace, disp,
this->M_offset,
this->M_identity );
546 displayer->leaderPrint (
" Non-Linear S - updating non linear terms in the Jacobian Matrix (Distributed Holzapfel): \n");
548 for(
UInt i(0); i <
this->M_vectorInterpolated.size(); i++ )
551 displayer->leaderPrint (
" ", i + 1,
"-th fiber family \n" );
554 Real kappaI =
this->M_dataMaterial->ithDistributionFibers( i );
555 Real alphaI =
this->M_dataMaterial->ithStiffnessFibers( i );
556 Real gammaI =
this->M_dataMaterial->ithNonlinearityFibers( i );
560 interpolatedValue_Type fiberIth = ExpressionDefinitions::interpolateFiber(
this->M_dispETFESpace, *(
this->M_vectorInterpolated[ i ] ) );
575 distributedStretch_Type distrStretch = ExpressionDistributedModel::distributedStretch( trCbar, IVithBar, kappaI );
578 tensorialPart_distrType tensorialPart = ExpressionDistributedModel::tensorialPartPiola( kappaI, I_C, IVith, F, F_T, Mith );
581 derivativeDistrStretch_Type derivativeStretch = ExpressionDistributedModel::derivativeDistributedStretch( kappaI, trCbar, IVithBar,
585 integrate( elements (
this->M_dispETFESpace->mesh() ),
586 this->M_dispFESpace->qr(),
587 this->M_dispETFESpace,
588 this->M_dispETFESpace,
589 value( 2.0 * alphaI ) * exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) *
590 Jel * distrStretch * derAtan( distrStretch ,
this->M_epsilon, ( 1.0 /
PI ) ) *
591 derivativeStretch * dot( tensorialPart , grad(phi_i) )
595 integrate( elements (
this->M_dispETFESpace->mesh() ),
596 this->M_dispFESpace->qr(),
597 this->M_dispETFESpace,
598 this->M_dispETFESpace,
599 atan( distrStretch ,
this->M_epsilon, ( 1 /
PI ), ( 1.0/2.0 ) ) *
600 value( 2.0 * alphaI ) * exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) *
601 Jel * derivativeStretch * dot( tensorialPart , grad(phi_i) )
605 integrate( elements (
this->M_dispETFESpace->mesh() ),
606 this->M_dispFESpace->qr(),
607 this->M_dispETFESpace,
608 this->M_dispETFESpace,
609 atan( distrStretch ,
this->M_epsilon, ( 1 /
PI ), ( 1.0/2.0 ) ) *
610 value( 2.0 * alphaI ) * Jel * distrStretch *
611 exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) *
612 value( 2.0 * gammaI ) * distrStretch * derivativeStretch * dot( tensorialPart , grad(phi_i) )
616 integrate( elements (
this->M_dispETFESpace->mesh() ),
617 this->M_dispFESpace->qr(),
618 this->M_dispETFESpace,
619 this->M_dispETFESpace,
620 atan( distrStretch ,
this->M_epsilon, ( 1 /
PI ), ( 1.0/2.0 ) ) *
621 value( 2.0 * alphaI ) * distrStretch *
622 exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) *
623 value( -2.0/3.0 ) * Jel * dot( F_T, grad(phi_j) ) * dot( tensorialPart , grad(phi_i) )
627 integrate( elements (
this->M_dispETFESpace->mesh() ),
628 this->M_dispFESpace->qr(),
629 this->M_dispETFESpace,
630 this->M_dispETFESpace,
631 atan( distrStretch ,
this->M_epsilon, ( 1 /
PI ), ( 1.0/2.0 ) ) *
632 value( 2.0 * alphaI ) * distrStretch *
633 exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) * Jel *
634 dot( value( kappaI ) * grad(phi_j) +
635 value( - ( 2.0 * kappaI )/ 3.0 ) * dot( F, grad(phi_j) ) * F_T +
636 value( kappaI / 3.0 ) * I_C * F_T * transpose( grad(phi_j) ) * F_T +
637 value( 1.0 - 3.0 * kappaI ) * grad(phi_j) * Mith +
638 value( ( 1.0 - 3.0*kappaI ) / 3.0 ) * IVith * F_T * transpose( grad(phi_j) ) * F_T +
639 value(-( 1.0 - 3.0*kappaI ) / 3.0 ) * ( dot( transpose(grad(phi_j)) * F , Mith ) + dot( transpose(F)*grad(phi_j), Mith ) ) * F_T , grad(phi_i) )
646 jacobian->globalAssemble();
653 template <
typename MeshType>
668 displayer->leaderPrint (
" \n*********************************\n ");
669 displayer->leaderPrint (
" Non-Linear S- Computing the Distributed Holzapfel nonlinear stiffness vector (Distributed Holzapfel) ");
670 displayer->leaderPrint (
" \n*********************************\n ");
681 tensorF_Type F = ExpressionDefinitions::deformationGradient(
this->M_dispETFESpace, disp,
this->M_offset,
this->M_identity );
701 for(
UInt i(0); i <
this->M_vectorInterpolated.size() ; i++ )
708 interpolatedValue_Type fiberIth = ExpressionDefinitions::interpolateFiber(
this->M_dispETFESpace, *(
this->M_vectorInterpolated[ i ] ) );
722 Real kappaI =
this->M_dataMaterial->ithDistributionFibers( i );
723 Real alphaI =
this->M_dataMaterial->ithStiffnessFibers( i );
724 Real gammaI =
this->M_dataMaterial->ithNonlinearityFibers( i );
727 distributedStretch_Type distrStretch = ExpressionDistributedModel::distributedStretch( trCbar, IVithBar, kappaI );
743 tensorialPart_distrType tensorialPart = ExpressionDistributedModel::tensorialPartPiola( kappaI, I_C, IVith, F, F_T, Mith );
745 integrate ( elements (
this->M_dispETFESpace->mesh() ),
746 this->M_dispFESpace->qr(),
747 this->M_dispETFESpace,
748 atan( distrStretch ,
this->M_epsilon, ( 1 /
PI ), ( 1.0/2.0 ) ) *
749 value( 2.0 * alphaI ) * Jel * ( distrStretch ) *
750 exp( value( gammaI ) * ( distrStretch ) * ( distrStretch ) ) *
751 dot( tensorialPart, grad( phi_i ) )
756 this->M_stiff->globalAssemble();
760 template <
typename MeshType>
769 template <
typename MeshType>
780 template <
typename MeshType>
782 const Epetra_SerialDenseMatrix& tensorF,
783 const Epetra_SerialDenseMatrix& cofactorF,
784 const std::vector<Real>& invariants,
791 template <
typename MeshType>
799 ASSERT( 2 < 0 ,
"For the distributed anisotropic law the computation of the Cauchy stress has to be defined." );
804 template <
typename MeshType>
super::isochoricStretch_Type isochoricStretch_Type
StructuralAnisotropicConstitutiveLaw< MeshType > * createDistributedHolzapfelMaterialNonLinear()
super::matrix_Type matrix_Type
void assignFunction(bcBase_Type &base)
Assign the function to the base of the BCHandler.
matrixSmall_Type M_identity
super::traceTensor_Type traceTensor_Type
ExpressionAddition< ExpressionAddition< scaledTensorF_Type, scaledTraceTimesMinusTF_Type >, ExpressionAddition< scaledFtimesM_Type, scaledFourthInvariantTimesMinusTF_Type > > tensorialPart_distrType
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.
void computeKinematicsVariables(const VectorElemental &dk_loc)
Computes the new Stiffness vector for Neo-Hookean and Holzapfel materials in StructuralSolver given a...
super::vectorFiberFunction_Type vectorFiberFunction_Type
super::vectorFiberFunctionPtr_Type vectorFiberFunctionPtr_Type
virtual ~DistributedHolzapfelMaterialNonLinear()
super::isochoricTrace_Type isochoricTrace_Type
vectorPtr_Type const activationDisplacement(const UInt) const
Get the activation displacement.
super::matrixPtr_Type matrixPtr_Type
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.
matrixPtr_Type const stiffMatrix() const
Get the Stiffness matrix.
DistributedHolzapfelMaterialNonLinear()
super::vectorsParametersPtr_Type vectorsParametersPtr_Type
vectorPtr_Type const stiffVector() const
Get the stiffness vector.
super::ETFESpacePtr_Type ETFESpacePtr_Type
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)...
super::FESpacePtr_Type FESpacePtr_Type
static const LifeV::UInt elm_nodes_num[]
super::vector_Type vector_Type
super::fiberFunction_Type fiberFunction_Type
mapMarkerIndexes_Type::const_iterator mapIteratorIndex_Type
super::data_Type data_Type
super::mapMarkerIndexesPtr_Type mapMarkerIndexesPtr_Type
void evaluateActivationFibers(const vector_Type &displacement, vector_Type &fourthInvariant)
ExpressionDistributedModel::linearizationDistributedStretch_Type derivativeDistrStretch_Type
super::vectorPtr_Type vectorPtr_Type
super::vectorsParameters_Type vectorsParameters_Type
void computeLinearStiff(dataPtr_Type &, const mapMarkerVolumesPtr_Type, const mapMarkerIndexesPtr_Type)
Compute the Stiffness matrix in StructuralSolver::buildSystem()
void setupFiberDirections(vectorFiberFunctionPtr_Type vectorOfFibers)
std::shared_ptr< vectorVolumes_Type > vectorVolumesPtr_Type
vectorPtr_Type const activatedDeterminant(const UInt) const
Get the activation displacement.
super::tensorF_Type tensorF_Type
super::FESpace_Type FESpace_Type
StructuralAnisotropicConstitutiveLaw< MeshType > super
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.
std::vector< UInt > vectorIndexes_Type
std::vector< typename MeshType::element_Type * > vectorVolumes_Type
void apply(const vector_Type &sol, vector_Type &res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const mapMarkerIndexesPtr_Type mapsMarkerIndexes, const displayerPtr_Type displayer)
ExpressionAddition< ExpressionProduct< ExpressionScalar, linearizationFisochoricTrace_Type >, ExpressionProduct< ExpressionScalar, linearizationFisochoricFourthInvariant_Type > > linearizationDistributedStretch_Type
super::displayerPtr_Type displayerPtr_Type
super::mapMarkerIndexes_Type mapMarkerIndexes_Type
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.
super::stretch_Type stretch_Type
vectorPtr_Type const selectionCriterion(const UInt) const
Specific for multimechanism.
super::powerExpression_Type powerExpression_Type
double Real
Generic real data.
super::determinantF_Type determinantF_Type
End namespace ExpressionDefinitions.
void computeReferenceConfigurations(const vector_Type &disp, const dataPtr_Type &dataMaterial, const displayerPtr_Type &displayer)
super::outerProduct_Type outerProduct_Type
super::minusT_Type minusT_Type
vectorPtr_Type const activatedUnitFiber(const UInt) const
Get the activation displacement.
std::shared_ptr< vectorIndexes_Type > vectorIndexesPtr_Type
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.
super::mapMarkerVolumes_Type mapMarkerVolumes_Type
QuadratureRule - The basis class for storing and accessing quadrature rules.
super::interpolatedValue_Type interpolatedValue_Type
vectorPtr_Type M_stiff
construct the vectors for the parameters
super::tensorC_Type tensorC_Type
super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type
uint32_type UInt
generic unsigned integer (used mainly for addressing)
mapMarkerVolumes_Type::const_iterator mapIterator_Type
super::dataPtr_Type dataPtr_Type
ExpressionAddition< distributedInvariants_Type, ExpressionScalar > distributedStretch_Type
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...
super::fiberFunctionPtr_Type fiberFunctionPtr_Type
MatrixSmall< 3, 3 > matrixSmall_Type