Commit 37b4fa0c authored by Paolo Tricerri's avatar Paolo Tricerri
Browse files

splitted long lines in the material class

parent 01606d3c
......@@ -67,16 +67,6 @@ public:
typedef typename super::dataPtr_Type dataPtr_Type;
typedef typename super::displayerPtr_Type displayerPtr_Type;
//! Definition of local tensors
typedef KN<Real> KN_Type;
typedef boost::shared_ptr<KN_Type> KNPtr_Type;
typedef KNM<Real> KNM_Type;
typedef boost::shared_ptr<KNM_Type> KNMPtr_Type;
typedef KNMK<Real> KNMK_Type;
typedef boost::shared_ptr<KNMK_Type> KNMKPtr_Type;
typedef typename super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type;
typedef typename super::mapMarkerVolumes_Type mapMarkerVolumes_Type;
typedef typename mapMarkerVolumes_Type::const_iterator mapIterator_Type;
......@@ -119,7 +109,8 @@ public:
//! Updates the Jacobian matrix in StructualSolver::updateJacobian
/*!
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateJacobianMatrix( const vector_Type& disp,
......@@ -132,7 +123,8 @@ public:
/*!
\param stiff: stiffness matrix provided from outside
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateNonLinearJacobianTerms( matrixPtr_Type& jacobian,
......@@ -147,18 +139,22 @@ public:
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the
material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const displayerPtr_Type& displayer );
//! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in StructuralSolver given a certain displacement field.
//! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
//! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in StructuralSolver
//! given a certain displacement field.
//! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem
//! since the matrix is the expression of the matrix is the same.
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void computeVector( const vector_Type& sol,
......@@ -168,21 +164,14 @@ public:
const displayerPtr_Type& displayer );
//! Computes the deformation gradient F, the cofactor matrix Cof(F), the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
//! Computes the deformation gradient F, the cofactor matrix Cof(F),
//! the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
//! This function is used in StructuralConstitutiveLaw::computeStiffness
/*!
\param dk_loc: the elemental displacement
*/
void computeKinematicsVariables( const VectorElemental& dk_loc );
//! Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F), the trace of C Tr(C).
/*!
\param dk_loc: local displacement vector
*/
//void computeStress( const vector_Type& sol );
//! ShowMe method of the class (saved on a file the stiffness vector and the jacobian)
void showMe( std::string const& fileNameVectStiff,
std::string const& fileNameJacobain );
......@@ -199,7 +188,8 @@ public:
//! Get the stiffness vector
vectorPtr_Type const stiffVector() const {return M_stiff; }
void apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes) ;
void apply( const vector_Type& sol, vector_Type& res,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes) ;
//@}
......@@ -259,7 +249,8 @@ template <typename Mesh>
void
ExponentialMaterialNonLinear<Mesh>::setup( const boost::shared_ptr< FESpace<Mesh, MapEpetra> >& dFESpace,
const boost::shared_ptr<const MapEpetra>& monolithicMap,
const UInt offset, const dataPtr_Type& dataMaterial, const displayerPtr_Type& displayer )
const UInt offset,
const dataPtr_Type& dataMaterial, const displayerPtr_Type& displayer )
{
this->M_displayer = displayer;
this->M_dataMaterial = dataMaterial;
......@@ -290,9 +281,10 @@ ExponentialMaterialNonLinear<Mesh>::setup( const boost::shared_ptr< FESpace<Mesh
template <typename Mesh>
void ExponentialMaterialNonLinear<Mesh>::computeLinearStiff(dataPtr_Type& /*dataMaterial*/, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/)
void ExponentialMaterialNonLinear<Mesh>::computeLinearStiff(dataPtr_Type& /*dataMaterial*/,
const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/)
{
//! Empty method for neo-hookean material
//! Empty method for exponential material
}
......@@ -381,41 +373,54 @@ void ExponentialMaterialNonLinear<Mesh>::updateNonLinearJacobianTerms( matrixPtr
//! VOLUMETRIC PART
//! 1. Stiffness matrix: int { 1/2 * bulk * ( 2 - 1/J + 1/J^2 ) * ( CofF : \nabla \delta ) (CofF : \nabla v) }
AssemblyElementalStructure::stiff_Jac_Pvol_1term( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_Pvol_1term( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! 2. Stiffness matrix: int { 1/2 * bulk * ( 1/J- 1 - log(J)/J^2 ) * ( CofF1 [\nabla \delta]^t CofF ) : \nabla v }
AssemblyElementalStructure::stiff_Jac_Pvol_2term( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_Pvol_2term( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! ISOCHORIC PART
//! 1. Stiffness matrix : int { - 2/3 * alpha * J^(-5/3) * exp( gamma*( Ic_iso - 3) )* ( 1. + coefExp * Ic_iso )
//! *( CofF : \nabla \delta ) ( F : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_1term( (-2.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_1term( (-2.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 2. Stiffness matrix : int { 2 * alpha * gamma * J^(-4/3) * exp( gamma*( Ic_iso - 3) ) *
//! ( F : \nabla \delta ) ( F : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_2term( 2.0 * alpha * gamma, gamma, (*M_Fk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_2term( 2.0 * alpha * gamma, gamma, (*M_Fk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 3. Stiffness matrix : int { 2.0/9.0 * alpha * J^-2 * Ic_iso * exp( gamma*( Ic_iso - 3) )*
//! ( 1. + gamma * Ic_iso )( CofF : \nabla \delta ) ( CofF : \nabla \v )}
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_3term( (2.0/9.0) * alpha, gamma, (*M_CofFk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_3term( (2.0/9.0) * alpha, gamma, (*M_CofFk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 4. Stiffness matrix : int { -2.0/3.0 * alpha * J^(-5/3) * exp( gamma*( Ic_iso - 3) )
//! * ( 1. + gamma * Ic_iso )( F : \nabla \delta ) ( CofF : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_4term( (-2.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_4term( (-2.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 5. Stiffness matrix : int { alpha * J^(-2/3) * exp( gamma*( Ic_iso - 3)) (\nabla \delta: \nabla \v)}
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_5term( alpha, gamma, (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_5term( alpha, gamma, (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 6. Stiffness matrix : int { 1.0/3.0 * alpha * J^(-2) * Ic_iso * exp(gamma*( Ic_iso - 3)) *
//! (CofF [\nabla \delta]^t CofF ) : \nabla \v }
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_6term( (1.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_Exp_6term( (1.0/3.0) * alpha, gamma, (*M_CofFk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! assembling
for ( UInt ic = 0; ic < nc; ++ic )
{
for ( UInt jc = 0; jc < nc; jc++ )
{
assembleMatrix( *jacobian, *this->M_elmatK, this->M_FESpace->fe(), this->M_FESpace->dof(), ic, jc, this->M_offset + ic*totalDof, this->M_offset + jc*totalDof );
assembleMatrix( *jacobian,
*this->M_elmatK,
this->M_FESpace->fe(),
this->M_FESpace->dof(),
ic, jc,
this->M_offset + ic*totalDof, this->M_offset + jc*totalDof );
}
}
}
......@@ -488,14 +493,16 @@ void ExponentialMaterialNonLinear<Mesh>::computeStiffness( const vector_Type& so
/*!
Source term Pvol: int { bulk /2* (J1^2 - J1 + log(J1) ) * 1/J1 * (CofF1 : \nabla v) }
*/
AssemblyElementalStructure::source_Pvol( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elvecK, this->M_FESpace->fe() );
AssemblyElementalStructure::source_Pvol( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elvecK, this->M_FESpace->fe() );
//! Isochoric part
/*!
Source term P1iso_Exp: int { alpha * exp(gamma *( Ic1_iso -3 )) *
( J1^(-2/3)* (F1 : \nabla v) - 1/3 * (Ic1_iso / J1) * (CofF1 : \nabla v) ) }
*/
AssemblyElementalStructure::source_P1iso_Exp( alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok), *this->M_elvecK, this->M_FESpace->fe() );
AssemblyElementalStructure::source_P1iso_Exp( alpha, gamma, (*M_CofFk), (*M_Fk), (*M_Jack), (*M_trCisok),
*this->M_elvecK, this->M_FESpace->fe() );
for ( UInt ic = 0; ic < nDimensions; ++ic )
{
......@@ -503,7 +510,10 @@ void ExponentialMaterialNonLinear<Mesh>::computeStiffness( const vector_Type& so
M_elvecK is assemble into *vec_stiff vector that is recall
from updateSystem(matrix_ptrtype& mat_stiff, vector_ptr_type& vec_stiff)
*/
assembleVector( *this->M_stiff, *this->M_elvecK, this->M_FESpace->fe(), this->M_FESpace->dof(), ic, this->M_offset + ic*totalDof );
assembleVector( *this->M_stiff,
*this->M_elvecK,
this->M_FESpace->fe(),
this->M_FESpace->dof(), ic, this->M_offset + ic*totalDof );
}
}
}
......@@ -534,7 +544,8 @@ void ExponentialMaterialNonLinear<Mesh>::computeKinematicsVariables( const Vecto
for ( UInt i = 0; i < this->M_FESpace->fe().nbFEDof(); i++ )
{
//! \grad u^k at a quadrature point
s += this->M_FESpace->fe().phiDer( i, jcoor, ig ) * dk_loc[ i + icoor * this->M_FESpace->fe().nbFEDof() ];
s += this->M_FESpace->fe().phiDer( i, jcoor, ig ) *
dk_loc[ i + icoor * this->M_FESpace->fe().nbFEDof() ];
}
//! gradient of displacement
(*M_Fk)[ icoor ][ jcoor ][ig ] = s;
......@@ -619,7 +630,8 @@ void ExponentialMaterialNonLinear<Mesh>::showMe( std::string const& fileNameStif
template <typename Mesh>
void ExponentialMaterialNonLinear<Mesh>::apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes )
void ExponentialMaterialNonLinear<Mesh>::apply( const vector_Type& sol,
vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes )
{
computeStiffness(sol, 0, this->M_dataMaterial, mapsMarkerVolumes, this->M_displayer);
res += *M_stiff;
......
......@@ -68,16 +68,6 @@ public:
typedef typename super::dataPtr_Type dataPtr_Type;
typedef typename super::displayerPtr_Type displayerPtr_Type;
//! Definition of local tensors
typedef KN<Real> KN_Type;
typedef boost::shared_ptr<KN_Type> KNPtr_Type;
typedef KNM<Real> KNM_Type;
typedef boost::shared_ptr<KNM_Type> KNMPtr_Type;
typedef KNMK<Real> KNMK_Type;
typedef boost::shared_ptr<KNMK_Type> KNMKPtr_Type;
typedef typename super::mapMarkerVolumesPtr_Type mapMarkerVolumesPtr_Type;
typedef typename super::mapMarkerVolumes_Type mapMarkerVolumes_Type;
typedef typename mapMarkerVolumes_Type::const_iterator mapIterator_Type;
......@@ -120,7 +110,8 @@ public:
//! Updates the Jacobian matrix in StructualSolver::updateJacobian
/*!
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateJacobianMatrix( const vector_Type& disp,
......@@ -133,7 +124,8 @@ public:
/*!
\param stiff: stiffness matrix provided from outside
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateNonLinearJacobianTerms( matrixPtr_Type& jacobian,
......@@ -148,18 +140,24 @@ public:
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const displayerPtr_Type& displayer );
void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
const displayerPtr_Type& displayer );
//! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in StructuralSolver given a certain displacement field.
//! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
//! Computes the new Stiffness vector for Neo-Hookean and Exponential materials in
//! StructuralSolver given a certain displacement field.
//! This function is used both in StructuralSolver::evalResidual and in StructuralSolver::updateSystem
//! since the matrix is the expression of the matrix is the same.
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void computeVector( const vector_Type& sol,
......@@ -169,21 +167,14 @@ public:
const displayerPtr_Type& displayer );
//! Computes the deformation gradient F, the cofactor matrix Cof(F), the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
//! Computes the deformation gradient F, the cofactor matrix Cof(F),
//! the determinant of F (J = det(F)), the trace of right Cauchy-Green tensor tr(C)
//! This function is used in StructuralConstitutiveLaw::computeStiffness
/*!
\param dk_loc: the elemental displacement
*/
void computeKinematicsVariables( const VectorElemental& dk_loc );
//! Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F), the trace of C Tr(C).
/*!
\param dk_loc: local displacement vector
*/
//void computeStress( const vector_Type& sol);
//! ShowMe method of the class (saved on a file the stiffness vector and the jacobian)
void showMe( std::string const& fileNameVectStiff,
std::string const& fileNameJacobain);
......@@ -201,7 +192,8 @@ public:
//! Get the stiffness vector
vectorPtr_Type const stiffVector() const {return M_stiff; }
void apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes);
void apply( const vector_Type& sol, vector_Type& res,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes);
//@}
......@@ -297,7 +289,8 @@ NeoHookeanMaterialNonLinear<Mesh>::setup( const boost::shared_ptr< FESpace<Mesh,
template <typename Mesh>
void NeoHookeanMaterialNonLinear<Mesh>::computeLinearStiff(dataPtr_Type& /*dataMaterial*/, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/)
void NeoHookeanMaterialNonLinear<Mesh>::computeLinearStiff(dataPtr_Type& /*dataMaterial*/,
const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/)
{
//! Empty method for neo-hookean material
}
......@@ -386,33 +379,44 @@ void NeoHookeanMaterialNonLinear<Mesh>::updateNonLinearJacobianTerms( matrixPtr_
//! VOLUMETRIC PART
//! 1. Stiffness matrix: int { 1/2 * bulk * ( 2 - 1/J + 1/J^2 ) * ( CofF : \nabla \delta ) (CofF : \nabla v) }
AssemblyElementalStructure::stiff_Jac_Pvol_1term( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_Pvol_1term( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! 2. Stiffness matrix: int { 1/2 * bulk * ( 1/J- 1 - log(J)/J^2 ) * ( CofF [\nabla \delta]^t CofF ) : \nabla v }
AssemblyElementalStructure::stiff_Jac_Pvol_2term( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_Pvol_2term( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! ISOCHORIC PART
//! 1. Stiffness matrix : int { -2/3 * mu * J^(-5/3) *( CofF : \nabla \delta ) ( F : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_NH_1term( (-2.0/3.0) * mu, (*M_CofFk), (*M_Fk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_NH_1term( (-2.0/3.0) * mu, (*M_CofFk), (*M_Fk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! 2. Stiffness matrix : int { 2/9 * mu * ( Ic_iso / J^2 )( CofF : \nabla \delta ) ( CofF : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_NH_2term( (2.0/9.0) * mu, (*M_CofFk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_NH_2term( (2.0/9.0) * mu, (*M_CofFk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! 3. Stiffness matrix : int { mu * J^(-2/3) (\nabla \delta : \nabla \v)}
AssemblyElementalStructure::stiff_Jac_P1iso_NH_3term( mu, (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
//! 4. Stiffness matrix : int { -2/3 * mu * J^(-5/3) ( F : \nabla \delta ) ( CofF : \nabla \v ) }
AssemblyElementalStructure::stiff_Jac_P1iso_NH_4term( (-2.0/3.0) * mu, (*M_CofFk), (*M_Fk), (*M_Jack), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_NH_4term( (-2.0/3.0) * mu, (*M_CofFk), (*M_Fk), (*M_Jack),
*this->M_elmatK, this->M_FESpace->fe() );
//! 5. Stiffness matrix : int { 1/3 * mu * J^(-2) * Ic_iso * (CofF [\nabla \delta]^t CofF ) : \nabla \v }
AssemblyElementalStructure::stiff_Jac_P1iso_NH_5term( (1.0/3.0) * mu, (*M_CofFk), (*M_Jack), (*M_trCisok), *this->M_elmatK, this->M_FESpace->fe() );
AssemblyElementalStructure::stiff_Jac_P1iso_NH_5term( (1.0/3.0) * mu, (*M_CofFk), (*M_Jack), (*M_trCisok),
*this->M_elmatK, this->M_FESpace->fe() );
//! assembling
for ( UInt ic = 0; ic < nc; ++ic )
{
for ( UInt jc = 0; jc < nc; jc++ )
{
assembleMatrix( *jacobian, *this->M_elmatK, this->M_FESpace->fe(), this->M_FESpace->dof(), ic, jc, this->M_offset + ic*totalDof, this->M_offset + jc*totalDof );
assembleMatrix( *jacobian,
*this->M_elmatK,
this->M_FESpace->fe(),
this->M_FESpace->dof(),
ic, jc,
this->M_offset + ic*totalDof, this->M_offset + jc*totalDof );
}
}
}
......@@ -420,7 +424,8 @@ void NeoHookeanMaterialNonLinear<Mesh>::updateNonLinearJacobianTerms( matrixPtr_
}
template <typename Mesh>
void NeoHookeanMaterialNonLinear<Mesh>::apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes )
void NeoHookeanMaterialNonLinear<Mesh>::apply( const vector_Type& sol, vector_Type& res,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes )
{
computeStiffness(sol, 0., this->M_dataMaterial, mapsMarkerVolumes, this->M_displayer);
res += *M_stiff;
......@@ -488,13 +493,15 @@ void NeoHookeanMaterialNonLinear<Mesh>::computeStiffness( const vector_Type&
/*!
Source term Pvol: int { bulk /2* (J1^2 - J1 + log(J1) ) * 1/J1 * (CofF1 : \nabla v) }
*/
AssemblyElementalStructure::source_Pvol( 0.5 * bulk, (*M_CofFk), (*M_Jack), *this->M_elvecK, this->M_FESpace->fe() );
AssemblyElementalStructure::source_Pvol( 0.5 * bulk, (*M_CofFk), (*M_Jack),
*this->M_elvecK, this->M_FESpace->fe() );
//! Isochoric part
/*!
Source term P1iso_NH
*/
AssemblyElementalStructure::source_P1iso_NH( mu, (*M_CofFk) ,(*M_Fk), (*M_Jack), (*M_trCisok) , *this->M_elvecK, this->M_FESpace->fe());
AssemblyElementalStructure::source_P1iso_NH( mu, (*M_CofFk) ,(*M_Fk), (*M_Jack), (*M_trCisok) ,
*this->M_elvecK, this->M_FESpace->fe());
for ( UInt ic = 0; ic < nDimensions; ++ic )
{
......@@ -502,7 +509,11 @@ void NeoHookeanMaterialNonLinear<Mesh>::computeStiffness( const vector_Type&
M_elvecK is assemble into *vec_stiff vector that is recall
from updateSystem(matrix_ptrtype& mat_stiff, vector_ptr_type& vec_stiff)
*/
assembleVector( *this->M_stiff, *this->M_elvecK, this->M_FESpace->fe(), this->M_FESpace->dof(), ic, this->M_offset + ic*totalDof );
assembleVector( *this->M_stiff,
*this->M_elvecK,
this->M_FESpace->fe(),
this->M_FESpace->dof(),
ic, this->M_offset + ic*totalDof );
}
}
}
......@@ -532,7 +543,8 @@ void NeoHookeanMaterialNonLinear<Mesh>::computeKinematicsVariables( const Vector
for ( UInt i = 0; i < this->M_FESpace->fe().nbFEDof(); i++ )
{
//! \grad u^k at a quadrature point
s += this->M_FESpace->fe().phiDer( i, jcoor, ig ) * dk_loc[ i + icoor * this->M_FESpace->fe().nbFEDof() ];
s += this->M_FESpace->fe().phiDer( i, jcoor, ig ) *
dk_loc[ i + icoor * this->M_FESpace->fe().nbFEDof() ];
}
//! gradient of displacement
(*M_Fk)[ icoor ][ jcoor ][ig ] = s;
......
......@@ -137,36 +137,46 @@ public:
*/
virtual void setup( const boost::shared_ptr< FESpace<Mesh, MapEpetra> >& dFESpace,
const boost::shared_ptr<const MapEpetra>& monolithicMap,
const UInt offset, const dataPtr_Type& dataMaterial, const displayerPtr_Type& displayer )=0;
const UInt offset, const dataPtr_Type& dataMaterial,
const displayerPtr_Type& displayer )=0;
//! Computes the linear part of the stiffness matrix StructuralSolver::buildSystem
/*!
\param dataMaterial the class with Material properties data
*/
virtual void computeLinearStiff( dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes ) = 0;
virtual void computeLinearStiff( dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes ) = 0;
//! Updates the Jacobian matrix in StructuralSolver::updateJacobian
/*!
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the
material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
virtual void updateJacobianMatrix( const vector_Type& disp, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const displayerPtr_Type& displayer ) = 0;
virtual void updateJacobianMatrix( const vector_Type& disp, const dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
const displayerPtr_Type& displayer ) = 0;
//! Computes the new Stiffness matrix in StructuralSolver given a certain displacement field. This function is used both in StructuralSolver::evalResidual and in
//! Computes the new Stiffness matrix in StructuralSolver given a certain displacement field.
//! This function is used both in StructuralSolver::evalResidual and in
//! StructuralSolver::updateSystem since the matrix is the expression of the matrix is the same.
//!This is virtual and not pure virtual since in the linear St. Venant-Kirchhoff law it is not needed.
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the
material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
virtual void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const displayerPtr_Type& displayer ) = 0;
virtual void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
const displayerPtr_Type& displayer ) = 0;
//! Computes the deformation Gradient F, the cofactor of F Cof(F), the determinant of F J = det(F), the trace of C Tr(C).
//! Computes the deformation Gradient F, the cofactor of F Cof(F),
//! the determinant of F J = det(F), the trace of C Tr(C).
/*!
\param dk_loc: local displacement vector
*/
......@@ -208,7 +218,8 @@ public:
//! Get the Stiffness matrix
virtual vectorPtr_Type const stiffVector() const = 0;
virtual void apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type mapsMarkerVolumes) =0;
virtual void apply( const vector_Type& sol, vector_Type& res,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes) =0;
//@}
......
......@@ -102,7 +102,8 @@ public:
//! Updates the Jacobian matrix in StructualSolver::updateJacobian
/*!
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateJacobianMatrix( const vector_Type& disp,
......@@ -114,7 +115,8 @@ public:
/*!
\param stiff: stiffness matrix provided from outside
\param disp: solution at the k-th iteration of NonLinearRichardson Method
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get
the material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void updateNonLinearJacobianTerms( matrixPtr_Type& jacobian,
......@@ -128,13 +130,15 @@ public:
/*!
\param sol: the solution vector
\param factor: scaling factor used in FSI
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the material coefficients (e.g. Young modulus, Poisson ratio..)
\param dataMaterial: a pointer to the dataType member in StructuralSolver class to get the
material coefficients (e.g. Young modulus, Poisson ratio..)
\param displayer: a pointer to the Dysplaier member in the StructuralSolver class
*/
void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes, const displayerPtr_Type& displayer );
void computeStiffness( const vector_Type& sol, Real factor, const dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes,
const displayerPtr_Type& displayer );
//! Missing Documentation!!!
void computeKinematicsVariables( const VectorElemental& /*dk_loc*/ ){}
//! ShowMe method of the class (saved on a file the two matrices)
......@@ -155,7 +159,8 @@ public:
//! Get the Stiffness vector
vectorPtr_Type const stiffVector() const { vectorPtr_Type zero( new vector_Type()); return zero;}
void apply( const vector_Type& sol, vector_Type& res, const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/) {res += *M_stiff*sol;}
void apply( const vector_Type& sol, vector_Type& res,
const mapMarkerVolumesPtr_Type /*mapsMarkerVolumes*/) {res += *M_stiff*sol;}
//@}
......@@ -207,7 +212,8 @@ VenantKirchhoffMaterialLinear<Mesh>::setup(const boost::shared_ptr< FESpace<Mesh
}
template <typename Mesh>
void VenantKirchhoffMaterialLinear<Mesh>::computeLinearStiff(dataPtr_Type& dataMaterial, const mapMarkerVolumesPtr_Type mapsMarkerVolumes)
void VenantKirchhoffMaterialLinear<Mesh>::computeLinearStiff(dataPtr_Type& dataMaterial,
const mapMarkerVolumesPtr_Type mapsMarkerVolumes)
{
// std::cout<<"compute LinearStiff Matrix start\n";
......@@ -249,7 +255,14 @@ void VenantKirchhoffMaterialLinear<Mesh>::computeLinearStiff(dataPtr_Type& dataM
{
for ( UInt jc = 0; jc < nc; jc++ )
{
assembleMatrix( *this->M_linearStiff, *this->M_elmatK, this->M_FESpace->fe(), this->M_FESpace->fe(), this->M_FESpace->dof(), this->M_FESpace->dof(), ic, jc, this->M_offset +ic*totalDof, this->M_offset + jc*totalDof );
assembleMatrix( *this->M_linearStiff,
*this->M_elmatK,
this->M_FESpace->fe(),
this->M_FESpace->fe(),
this->M_FESpace->dof(),
this->M_FESpace->dof(),
ic, jc,
this->M_offset +ic*totalDof, this->M_offset + jc*totalDof );
}
}
......