Chaste  Release::3.4
AbstractNonlinearElasticitySolver< DIM > Class Template Referenceabstract

#include <AbstractNonlinearElasticitySolver.hpp>

+ Inheritance diagram for AbstractNonlinearElasticitySolver< DIM >:
+ Collaboration diagram for AbstractNonlinearElasticitySolver< DIM >:

Public Member Functions

void ComputeResidual (Vec currentGuess, Vec residualVector)
 
void ComputeJacobian (Vec currentGuess, Mat *pJacobian, Mat *pPreconditioner)
 
 AbstractNonlinearElasticitySolver (AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory, CompressibilityType compressibilityType)
 
virtual ~AbstractNonlinearElasticitySolver ()
 
void Solve (double tol=-1.0)
 
void SetIncludeActiveTension (bool includeActiveTension=true)
 
unsigned GetNumNewtonIterations ()
 
void SetWriteOutputEachNewtonIteration (bool writeOutputEachNewtonIteration=true)
 
void SetKspAbsoluteTolerance (double kspAbsoluteTolerance)
 
void SetTakeFullFirstNewtonStep (bool takeFullFirstStep=true)
 
void SetUsePetscDirectSolve (bool usePetscDirectSolve=true)
 
void SetCurrentTime (double time)
 
void CreateCmguiOutput ()
 
void WriteCurrentStrains (StrainType strainType, std::string fileName, int counterToAppend=-1)
 
void SetComputeAverageStressPerElementDuringSolve (bool setComputeAverageStressPerElement=true)
 
void WriteCurrentAverageElementStresses (std::string fileName, int counterToAppend=-1)
 
std::vector< c_vector< double,
DIM > > & 
rGetSpatialSolution ()
 
std::vector< c_vector< double,
DIM > > & 
rGetDeformedPosition ()
 
c_matrix< double, DIM, DIM > GetAverageStressPerElement (unsigned elementIndex)
 
- Public Member Functions inherited from AbstractContinuumMechanicsSolver< DIM >
 AbstractContinuumMechanicsSolver (AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, ContinuumMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory, CompressibilityType compressibilityType)
 
virtual ~AbstractContinuumMechanicsSolver ()
 
void WriteCurrentSpatialSolution (std::string fileName, std::string fileExtension, int counterToAppend=-1)
 
void WriteCurrentPressureSolution (int counterToAppend=-1)
 
void SetWriteOutput (bool writeOutput=true)
 
void CreateVtkOutput (std::string spatialSolutionName="Spatial solution")
 
std::vector< double > & rGetCurrentSolution ()
 
std::vector< double > & rGetPressures ()
 

Protected Member Functions

void AddStressToAverageStressPerElement (c_matrix< double, DIM, DIM > &rT, unsigned elementIndex)
 
virtual void SetKspSolverAndPcType (KSP solver)
 
virtual void AssembleSystem (bool assembleResidual, bool assembleLinearSystem)=0
 
virtual void FinishAssembleSystem (bool assembleResidual, bool assembleLinearSystem)
 
void GetElementCentroidStrain (StrainType strainType, Element< DIM, DIM > &rElement, c_matrix< double, DIM, DIM > &rDeformationGradient)
 
virtual void AddActiveStressAndStressDerivative (c_matrix< double, DIM, DIM > &rC, unsigned elementIndex, unsigned currentQuadPointGlobalIndex, c_matrix< double, DIM, DIM > &rT, FourthOrderTensor< DIM, DIM, DIM, DIM > &rDTdE, bool addToDTdE)
 
virtual void SetupChangeOfBasisMatrix (unsigned elementIndex, unsigned currentQuadPointGlobalIndex)
 
void AssembleOnBoundaryElement (BoundaryElement< DIM-1, DIM > &rBoundaryElement, c_matrix< double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE > &rAelem, c_vector< double, BOUNDARY_STENCIL_SIZE > &rBelem, bool assembleResidual, bool assembleJacobian, unsigned boundaryConditionIndex)
 
bool ShouldAssembleMatrixTermForPressureOnDeformedBc ()
 
void AssembleOnBoundaryElementForPressureOnDeformedBc (BoundaryElement< DIM-1, DIM > &rBoundaryElement, c_matrix< double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE > &rAelem, c_vector< double, BOUNDARY_STENCIL_SIZE > &rBelem, bool assembleResidual, bool assembleJacobian, unsigned boundaryConditionIndex)
 
double ComputeResidualAndGetNorm (bool allowException)
 
double CalculateResidualNorm ()
 
void VectorSum (std::vector< double > &rX, ReplicatableVector &rY, double a, std::vector< double > &rZ)
 
void PrintLineSearchResult (double s, double residNorm)
 
double TakeNewtonStep ()
 
double UpdateSolutionUsingLineSearch (Vec solution)
 
virtual void PostNewtonStep (unsigned counter, double normResidual)
 
void SolveNonSnes (double tol=-1.0)
 
- Protected Member Functions inherited from AbstractContinuumMechanicsSolver< DIM >
void AllocateMatrixMemory ()
 
void ApplyDirichletBoundaryConditions (ApplyDirichletBcsType type, bool symmetricProblem)
 
void AddIdentityBlockForDummyPressureVariables (ApplyDirichletBcsType type)
 
void RemovePressureDummyValuesThroughLinearInterpolation ()
 

Protected Attributes

SolidMechanicsProblemDefinition
< DIM > & 
mrProblemDefinition
 
MatmrJacobianMatrix
 
c_matrix< double, DIM, DIM > mChangeOfBasisMatrix
 
double mKspAbsoluteTol
 
bool mWriteOutputEachNewtonIteration
 
FourthOrderTensor< DIM, DIM,
DIM, DIM > 
dTdE
 
unsigned mNumNewtonIterations
 
double mCurrentTime
 
bool mCheckedOutwardNormals
 
bool mUseSnesSolver
 
double mLastDampingValue
 
bool mFirstStep
 
bool mTakeFullFirstNewtonStep
 
bool mPetscDirectSolve
 
bool mIncludeActiveTension
 
bool mSetComputeAverageStressPerElement
 
std::vector< c_vector< double,
DIM *(DIM+1)/2 > > 
mAverageStressesPerElement
 
- Protected Attributes inherited from AbstractContinuumMechanicsSolver< DIM >
AbstractTetrahedralMesh< DIM,
DIM > & 
mrQuadMesh
 
ContinuumMechanicsProblemDefinition
< DIM > & 
mrProblemDefinition
 
bool mWriteOutput
 
std::string mOutputDirectory
 
OutputFileHandlermpOutputFileHandler
 
std::vector< c_vector< double,
DIM > > 
mSpatialSolution
 
std::vector< doublemPressureSolution
 
std::vector< doublemCurrentSolution
 
GaussianQuadratureRule< DIM > * mpQuadratureRule
 
GaussianQuadratureRule< DIM-1 > * mpBoundaryQuadratureRule
 
CompressibilityType mCompressibilityType
 
unsigned mProblemDimension
 
unsigned mNumDofs
 
bool mVerbose
 
Vec mResidualVector
 
Vec mLinearSystemRhsVector
 
Mat mSystemLhsMatrix
 
Vec mDirichletBoundaryConditionsVector
 
Mat mPreconditionMatrix
 

Static Protected Attributes

static const size_t NUM_VERTICES_PER_ELEMENT = DIM+1
 
static const size_t NUM_NODES_PER_ELEMENT = (DIM+1)*(DIM+2)/2
 
static const size_t NUM_NODES_PER_BOUNDARY_ELEMENT = DIM*(DIM+1)/2
 
static const size_t BOUNDARY_STENCIL_SIZE = DIM*NUM_NODES_PER_BOUNDARY_ELEMENT
 
static double MAX_NEWTON_ABS_TOL = 1e-7
 
static double MIN_NEWTON_ABS_TOL = 1e-10
 
static double NEWTON_REL_TOL = 1e-4
 

Private Member Functions

void SolveSnes ()
 

Friends

class StressRecoveror< DIM >
 
class VtkNonlinearElasticitySolutionWriter< DIM >
 

Detailed Description

template<unsigned DIM>
class AbstractNonlinearElasticitySolver< DIM >

Abstract nonlinear elasticity solver. IncompressibleNonlinearElasticityAssembler and CompressibleNonlinearElasticityAssembler inherit from this class.

The class is both a solver AND a assembler: the AssembleOnElement(), AssembleSystem() methods are hardcoded into this class. In principle something like AbstractContinuumMechanicsAssembler could have been used as a member variable [that class is used for assembling fluids matrices etc] but nonlinear elasticity is too complex for the that class to be used, as things like stress and stress-derivative need to be computed at the AssembleOnElement level, things like pressure-on-deformed-surface are deformation dependent boundary conditions, etc. *

Definition at line 139 of file AbstractNonlinearElasticitySolver.hpp.

Constructor & Destructor Documentation

template<unsigned DIM>
AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver ( AbstractTetrahedralMesh< DIM, DIM > &  rQuadMesh,
SolidMechanicsProblemDefinition< DIM > &  rProblemDefinition,
std::string  outputDirectory,
CompressibilityType  compressibilityType 
)

Constructor.

Parameters
rQuadMeshthe quadratic mesh
rProblemDefinitionan object defining in particular the body force and boundary conditions
outputDirectoryoutput directory
compressibilityTypeShould be equal to COMPRESSIBLE or INCOMPRESSIBLE (see enumeration defined at top of file) (depending on which concrete class is inheriting from this) and is only used in computing mNumDofs and allocating matrix memory.

Definition at line 823 of file AbstractNonlinearElasticitySolver.hpp.

References CommandLineArguments::Instance(), AbstractNonlinearElasticitySolver< DIM >::mChangeOfBasisMatrix, AbstractNonlinearElasticitySolver< DIM >::mPetscDirectSolve, AbstractNonlinearElasticitySolver< DIM >::mrProblemDefinition, AbstractNonlinearElasticitySolver< DIM >::mTakeFullFirstNewtonStep, AbstractNonlinearElasticitySolver< DIM >::mUseSnesSolver, and CommandLineArguments::OptionExists().

template<unsigned DIM>
AbstractNonlinearElasticitySolver< DIM >::~AbstractNonlinearElasticitySolver ( )
virtual

Destructor.

Definition at line 849 of file AbstractNonlinearElasticitySolver.hpp.

Member Function Documentation

template<unsigned DIM>
virtual void AbstractNonlinearElasticitySolver< DIM >::AddActiveStressAndStressDerivative ( c_matrix< double, DIM, DIM > &  rC,
unsigned  elementIndex,
unsigned  currentQuadPointGlobalIndex,
c_matrix< double, DIM, DIM > &  rT,
FourthOrderTensor< DIM, DIM, DIM, DIM > &  rDTdE,
bool  addToDTdE 
)
inlineprotectedvirtual

Add on the active component to the stress (and maybe also to the stress-derivative). This is called after getting the passive stress and stress-derivative from a material law.

This method does nothing but can be overloaded by the other solvers, see eg cardiac mechanics solvers.

Parameters
rCThe Lagrangian deformation tensor (F^T F)
elementIndexIndex of the current element
currentQuadPointGlobalIndexThe index (assuming an outer loop over elements and an inner loop over quadrature points), of the current quadrature point.
rTThe stress to be added to
rDTdEthe stress derivative to be added to, assuming the final parameter is true
addToDTdEA boolean flag saying whether the stress derivative is required or not.

Definition at line 429 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::AddStressToAverageStressPerElement ( c_matrix< double, DIM, DIM > &  rT,
unsigned  elementIndex 
)
protected

Add the given stress tensor to the store of average stresses. mSetComputeAverageStressPerElement must be true

Parameters
rT2nd PK stress (matrix is assumed symmetric)
elementIndexelement index

Definition at line 1008 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::AssembleOnBoundaryElement ( BoundaryElement< DIM-1, DIM > &  rBoundaryElement,
c_matrix< double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE > &  rAelem,
c_vector< double, BOUNDARY_STENCIL_SIZE > &  rBelem,
bool  assembleResidual,
bool  assembleJacobian,
unsigned  boundaryConditionIndex 
)
protected

Compute the term from the surface integral of s*phi, where s is a specified non-zero surface traction (ie Neumann boundary condition) to be added to the residual vector.

Calls AssembleOnBoundaryElementForPressureOnDeformedBc() if appropriate

Parameters
rBoundaryElementthe boundary element to be integrated on
rAelemThe element's contribution to the LHS matrix is returned. There is no need to zero this matrix before calling.
rBelemThe element's contribution to the RHS vector is returned. There is no need to zero this vector before calling.
assembleResidualA bool stating whether to assemble the residual vector.
assembleJacobianA bool stating whether to assemble the Jacobian matrix.
boundaryConditionIndexindex of this boundary (in the vectors in the problem definition object, in which the boundary conditions are stored

Definition at line 1185 of file AbstractNonlinearElasticitySolver.hpp.

References QuadraticBasisFunction< ELEMENT_DIM >::ComputeBasisFunctions(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetIndex(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetNodeGlobalIndex(), and NEVER_REACHED.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::AssembleOnBoundaryElementForPressureOnDeformedBc ( BoundaryElement< DIM-1, DIM > &  rBoundaryElement,
c_matrix< double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE > &  rAelem,
c_vector< double, BOUNDARY_STENCIL_SIZE > &  rBelem,
bool  assembleResidual,
bool  assembleJacobian,
unsigned  boundaryConditionIndex 
)
protected

Alternative version of AssembleOnBoundaryElement which is used for problems where a normal pressure is applied to the deformed body. The traction then depends on the current deformation, specifically s = -det(F)*P*F^{-T}N.

To compute F we have to find the volume element containing this surface element and use this in the computation. See comments in implementation for more details.

Parameters
rBoundaryElementthe boundary element to be integrated on
rAelemThe element's contribution to the LHS matrix is returned. There is no need to zero this matrix before calling.
rBelemThe element's contribution to the RHS vector is returned. There is no need to zero this vector before calling.
assembleResidualA bool stating whether to assemble the residual vector.
assembleJacobianA bool stating whether to assemble the Jacobian matrix.
boundaryConditionIndexindex of this boundary (in the vectors in the problem definition object, in which the boundary conditions are stored

Definition at line 1287 of file AbstractNonlinearElasticitySolver.hpp.

References Element< ELEMENT_DIM, SPACE_DIM >::CalculateInterpolationWeights(), AbstractTetrahedralElement< ELEMENT_DIM, SPACE_DIM >::CalculateNormal(), QuadraticBasisFunction< ELEMENT_DIM >::ComputeBasisFunctions(), QuadraticBasisFunction< ELEMENT_DIM >::ComputeTransformedBasisFunctionDerivatives(), Determinant(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetIndex(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetNode(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetNodeGlobalIndex(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetNumNodes(), Inverse(), and NEVER_REACHED.

template<unsigned DIM>
virtual void AbstractNonlinearElasticitySolver< DIM >::AssembleSystem ( bool  assembleResidual,
bool  assembleLinearSystem 
)
protectedpure virtual

Assemble the residual vector and/or Jacobian matrix (using the current solution stored in mCurrentSolution, output going to mResidualVector and/or mrJacobianMatrix).

Must be overridden in concrete derived classes.

Parameters
assembleResidualA bool stating whether to assemble the residual vector.
assembleLinearSystemA bool stating whether to assemble the Jacobian matrix and the RHS vector of the linear system (which is based on the residual but could be slightly different due to the way dirichlet boundary conditions are applied to the linear system - see comments in ApplyDirichletBoundaryConditions).

Implemented in IncompressibleNonlinearElasticitySolver< DIM >, and CompressibleNonlinearElasticitySolver< DIM >.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::CalculateResidualNorm ( )
protected
Returns
|r|_2 / length(r), where r is the current residual vector.

Definition at line 1715 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::ComputeJacobian ( Vec  currentGuess,
Mat pJacobian,
Mat pPreconditioner 
)

Public method for computing the jacobian that will be called, effectively, by the SNES solver

Parameters
currentGuessInput, the current guess for the solution
pJacobianOutput, the jacobian matrix at this guess
pPreconditionerOutput, the preconditioner matrix

Definition at line 2279 of file AbstractNonlinearElasticitySolver.hpp.

References GenericEventHandler< 7, MechanicsEventHandler >::BeginEvent(), GenericEventHandler< 7, MechanicsEventHandler >::EndEvent(), and ReplicatableVector::GetSize().

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::ComputeResidual ( Vec  currentGuess,
Vec  residualVector 
)

Public method for computing the residual that will be called, effectively, by the SNES solver

Parameters
currentGuessInput, the current guess for the solution
residualVectorOutput, the residual vector given this guess.

Definition at line 2263 of file AbstractNonlinearElasticitySolver.hpp.

References ReplicatableVector::GetSize().

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::ComputeResidualAndGetNorm ( bool  allowException)
protected

Set up the residual vector (using the current solution), and get its scaled norm (Calculate |r|_2 / length(r), where r is residual vector).

Parameters
allowExceptionSometimes the current solution solution will be such that the residual vector cannot be computed, as (say) the material law will throw an exception as the strains are too large. If this bool is set to true, the exception will be caught, and DBL_MAX returned as the residual norm
Returns
residual norm

Definition at line 1689 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::CreateCmguiOutput ( )

Convert the output to Cmgui format (placed in a folder called cmgui in the output directory). Writes the original mesh as solution_0.exnode and the (current) solution as solution_1.exnode.

Definition at line 979 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION, CmguiDeformedSolutionsWriter< DIM >::WriteCmguiScript(), CmguiDeformedSolutionsWriter< DIM >::WriteDeformationPositions(), and CmguiDeformedSolutionsWriter< DIM >::WriteInitialMesh().

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::FinishAssembleSystem ( bool  assembleResidual,
bool  assembleLinearSystem 
)
protectedvirtual

To be called at the end of AssembleSystem. Calls (Petsc) assemble methods on the Vecs and Mat, and calls ApplyDirichletBoundaryConditions.

Parameters
assembleResidualsee documentation for AssembleSystem
assembleLinearSystemsee documentation for AssembleSystem

Definition at line 855 of file AbstractNonlinearElasticitySolver.hpp.

References PetscVecTools::Finalise(), PetscMatTools::Finalise(), and PetscMatTools::SwitchWriteMode().

template<unsigned DIM>
c_matrix< double, DIM, DIM > AbstractNonlinearElasticitySolver< DIM >::GetAverageStressPerElement ( unsigned  elementIndex)

If SetComputeAverageStressPerElementDuringSolve() was called before the Solve(), then this method can be used to get the average stress for a particular element.

Parameters
elementIndexelementIndex
Returns
stress tensor

Definition at line 1042 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::GetElementCentroidStrain ( StrainType  strainType,
Element< DIM, DIM > &  rElement,
c_matrix< double, DIM, DIM > &  rDeformationGradient 
)
protected

Compute the strain at the centroid at an element

Parameters
strainTypeWhich strain to compute, should be one of: DEFORMATION_GRADIENT_F, DEFORMATION_TENSOR_C, or LAGRANGE_STRAIN_E
rElementThe element
rDeformationGradientReference to a matrix, which will be filled in by this method.

Definition at line 1085 of file AbstractNonlinearElasticitySolver.hpp.

References QuadraticBasisFunction< ELEMENT_DIM >::ComputeTransformedBasisFunctionDerivatives(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetIndex(), AbstractElement< ELEMENT_DIM, SPACE_DIM >::GetNodeGlobalIndex(), NEVER_REACHED, and ChastePoint< DIM >::rGetLocation().

template<unsigned DIM>
unsigned AbstractNonlinearElasticitySolver< DIM >::GetNumNewtonIterations ( )
Returns
number of Newton iterations taken in last solve.

Definition at line 2161 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::PostNewtonStep ( unsigned  counter,
double  normResidual 
)
protectedvirtual

This function may be overloaded by subclasses. It is called after each Newton iteration.

Parameters
countercurrent newton iteration number
normResidualnorm of the residual

Definition at line 2068 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::PrintLineSearchResult ( double  s,
double  residNorm 
)
protected

Print to std::cout the residual norm for this s, ie ||f(x+su)|| where f is the residual vector, x the current solution and u the update vector.

Parameters
ss
residNormresidual norm.

Definition at line 1918 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
std::vector< c_vector< double, DIM > > & AbstractNonlinearElasticitySolver< DIM >::rGetDeformedPosition ( )
Returns
the deformed position. Note: return_value[i](j) = x_j for node i. Just calls rGetSpatialSolution().

Definition at line 906 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
std::vector< c_vector< double, DIM > > & AbstractNonlinearElasticitySolver< DIM >::rGetSpatialSolution ( )
virtual
Returns
the deformed position. Note: return_value[i](j) = x_j for node i.

Implements AbstractContinuumMechanicsSolver< DIM >.

Definition at line 891 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetComputeAverageStressPerElementDuringSolve ( bool  setComputeAverageStressPerElement = true)

The user may request that the stress for each element (averaged over quadrature point stresses) are saved during the Solve(), by calling this.

Parameters
setComputeAverageStressPerElementwhether to compute stresses (defaults to true)

Definition at line 998 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetCurrentTime ( double  time)
inline

This solver is for static problems, however the body force or surface tractions could be a function of time. This method is for setting the time.

Parameters
timecurrent time

Definition at line 745 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mCurrentTime.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetIncludeActiveTension ( bool  includeActiveTension = true)
inline

Whether to call AddActiveStressAndStressDerivative() when computing stresses or not.

Subclasses, such as the cardiac mechanics solvers, may implement the above method to add active contributions to the stress. However, sometimes we might want to switch this off, which is what this function is for - will generally be called with includeActiveTension=false

Parameters
includeActiveTensionwhether to include active tension

Definition at line 667 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mIncludeActiveTension.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetKspAbsoluteTolerance ( double  kspAbsoluteTolerance)
inline

Set the absolute tolerance to be used when solving the linear system. If this is not called a relative tolerance is used.

Parameters
kspAbsoluteTolerancethe tolerance

Definition at line 696 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mKspAbsoluteTol.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetKspSolverAndPcType ( KSP  solver)
protectedvirtual

Set the KSP type (CG, GMRES, etc) and the preconditioner type (ILU, ICC etc). Depends on incompressible or not, and other factors.

///

Todo:
#2057 Make better choices here...
Parameters
solverKSP solver object (Petsc object)
Todo:
#2057 Make better choices..

Definition at line 1600 of file AbstractNonlinearElasticitySolver.hpp.

References PetscTools::IsSequential().

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetTakeFullFirstNewtonStep ( bool  takeFullFirstStep = true)
inline

The following odd behaviour has been observed: for some problems the solver will fail in the first Newton iteration, with the residual not decreasing in the direction of the Newton update, but: if you take a full Newton step anyway (increasing the residual-norm), the solver then converges perfectly. This method allows the user to choose this.

Does nothing if the SNES solver is used.

See ticket #2304

Parameters
takeFullFirstStepWhether to take a full first Newton step or not.

Definition at line 715 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mTakeFullFirstNewtonStep.

template<unsigned DIM>
virtual void AbstractNonlinearElasticitySolver< DIM >::SetupChangeOfBasisMatrix ( unsigned  elementIndex,
unsigned  currentQuadPointGlobalIndex 
)
inlineprotectedvirtual

Should re-set the variable mChangeOfBasisMatrix if this is not to be the identity. Here, does nothing, but over-ridden in cardiac mechanics solvers.

Parameters
elementIndexelement global index
currentQuadPointGlobalIndexglobal index of the quad point

Definition at line 446 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetUsePetscDirectSolve ( bool  usePetscDirectSolve = true)
inline

Get Petsc to do a direct solve on the linear system (instead of using an iterative solve). This is equivalent to passing in command line arguments -ksp_type pre_only -pc_type lu through to Petsc, but in the incompressible case the preconditioner is set equal to the Jacobian with a mass matrix in the pressure-pressure block (to avoid zeros on the diagonal. Hence a few linear solve iterations are required for this case. Using a direct solve can lead to huge computation time savings if there is enough memory for it: the linear solve may be faster and nonlinear convergence likely to be much better, as the linear solve is exact.

Parameters
usePetscDirectSolveWhether to use the Petsc direct solver or not

Definition at line 733 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mPetscDirectSolve.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetWriteOutputEachNewtonIteration ( bool  writeOutputEachNewtonIteration = true)
inline

By default only the original and converged solutions are written. Call this by get node positions written after every Newton step (mostly for debugging).

Parameters
writeOutputEachNewtonIterationwhether to write each iteration

Definition at line 685 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mWriteOutputEachNewtonIteration.

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::ShouldAssembleMatrixTermForPressureOnDeformedBc ( )
protected

When pressure-on-deformed-body boundary conditions are used:

For some reason the use of the Jacobian corresponding to the pressure term doesn't help (makes things worse!) if the current guess is not close enough to the solution, and can lead to immediate divergence. It will lead to faster convergence once close enough to the solution. This method contains the logic used to decide whether to include the jacobian for this term or not.

We only assemble the contribution to the matrix if the last damping value is close enough to 1 (non-snes solver). This will current always return false if the snes solver is being used

Returns
true if we are assembling pressure-on-deformed-body

Definition at line 1264 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::Solve ( double  tol = -1.0)

Solve the problem.

Parameters
toltolerance used in Newton solve (defaults to -1.0). Not used in SNES solves.

Definition at line 1553 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SolveNonSnes ( double  tol = -1.0)
protected

Solve method which uses a nonlinear solver coded in this class (as opposed to the SNES solver. Private, user should call Solve()

Parameters
tolabsolute solver used in nonlinear solve

Definition at line 2074 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SolveSnes ( )
private

Alternative solve method which uses a Petsc SNES solver. Private, user should call Solve()

Definition at line 2174 of file AbstractNonlinearElasticitySolver.hpp.

References PetscTools::Destroy(), EXCEPTION, PetscVecTools::Finalise(), and PETSC_DESTROY_PARAM.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::TakeNewtonStep ( )
protected

Take one newton step, by solving the linear system -Ju=f, (J the jacobian, f the residual, u the update), and picking s such that a_new = a_old + su (a the current solution) such |f(a)| is the smallest.

Returns
The current norm of the residual after the newton step.

Definition at line 1740 of file AbstractNonlinearElasticitySolver.hpp.

References GenericEventHandler< 7, MechanicsEventHandler >::BeginEvent(), PetscTools::Destroy(), GenericEventHandler< 7, MechanicsEventHandler >::EndEvent(), EXCEPTION, PetscTools::GetMyRank(), PETSC_DESTROY_PARAM, Timer::PrintAndReset(), and Timer::Reset().

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::UpdateSolutionUsingLineSearch ( Vec  solution)
protected

Using the update vector (of Newton's method), choose s such that ||f(x+su)|| is most decreased, where f is the residual vector, x the current solution (mCurrentSolution) and u the update vector. This checks s=1 first (most likely to be the current solution, then 0.9, 0.8.. until ||f|| starts increasing.

Parameters
solutionThe solution of the linear solve in newton's method, ie the update vector u.
Returns
norm of residual

Definition at line 1927 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::VectorSum ( std::vector< double > &  rX,
ReplicatableVector rY,
double  a,
std::vector< double > &  rZ 
)
protected

Simple helper function, computes Z = X + aY, where X and Z are std::vectors and Y a ReplicatableVector.

Parameters
rXX
rYY (replicatable vector)
aa
rZZ the returned vector

Definition at line 1725 of file AbstractNonlinearElasticitySolver.hpp.

References ReplicatableVector::GetSize().

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::WriteCurrentAverageElementStresses ( std::string  fileName,
int  counterToAppend = -1 
)

If SetComputeAverageStressPerElementDuringSolve() was called before the Solve(), then this method can be used to print the average stresses to file)

Each line of the output file corresponds to one element: the DIM*DIM matrix will be written as one line, using the ordering: T00 T01 T02 T10 T11 T12 T20 T21 T22.

Parameters
fileNameThe file name stem
counterToAppend(Optional) number to append in the filename.

The final file is [fileName]_[counterToAppend].stress

Definition at line 952 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::WriteCurrentStrains ( StrainType  strainType,
std::string  fileName,
int  counterToAppend = -1 
)

Write the strain for each element (evaluated at the centroids of each element). Which strain to compute is determined by the first input parameter, and will be either F, C or E. Each line of the output file corresponds to one element: the DIM*DIM matrix will be written as one line, using the following ordering (assuming F is written). F00 F01 F02 F10 F11 F12 F20 F21 F22.

Parameters
strainTypeWhich strain to write, should be one of: DEFORMATION_GRADIENT_F, DEFORMATION_TENSOR_C, or LAGRANGE_STRAIN_E
fileNameThe file name stem
counterToAppend(Optional) number to append in the filename.

The final file is [fileName]_[counterToAppend].strain

Definition at line 913 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM >::GetElementIteratorBegin().

Member Data Documentation

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::BOUNDARY_STENCIL_SIZE = DIM*NUM_NODES_PER_BOUNDARY_ELEMENT
staticprotected

Boundary stencil size. Note this is just the number of spatial unknowns on the boundary element, because the boundary integral terms (in either compressible or incompressible formulations) (i) do not involve pressure and (ii) do not appear in the pressure equations (the constraint equations).

Definition at line 223 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
FourthOrderTensor<DIM,DIM,DIM,DIM> AbstractNonlinearElasticitySolver< DIM >::dTdE
protected

Storage space for a 4th order tensor used in assembling the Jacobian (to avoid repeated memory allocation).

Definition at line 277 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
std::vector<c_vector<double,DIM*(DIM+1)/2> > AbstractNonlinearElasticitySolver< DIM >::mAverageStressesPerElement
protected

If the stresses for each element (averaged over quadrature point stresses) are to be stored, they are stored in this variable. Note to save memory we just don't store the lower half of the stress, as the stress is symmetric, hence this is a vector of 6 (in 3d) variables rather than a 3d matrix.

Definition at line 346 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::MAX_NEWTON_ABS_TOL = 1e-7
staticprotected

Maximum absolute tolerance for Newton solve. The Newton solver uses the absolute tolerance corresponding to the specified relative tolerance, but has a max and min allowable absolute tolerance. (ie: if max_abs = 1e-7, min_abs = 1e-10, rel=1e-4: then if the norm of the initial_residual (=a) is 1e-2, it will solve with tolerance 1e-7; if a=1e-5, it will solve with tolerance 1e-9; a=1e-9, it will solve with tolerance 1e-10.

Definition at line 232 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
c_matrix<double,DIM,DIM> AbstractNonlinearElasticitySolver< DIM >::mChangeOfBasisMatrix
protected

Matrix to be passed to material law, in case the material law is anisotropic and depends on a local coordinate system (eg cardiac problems). Defaults to the identity matrix. The cardiac mechanics solvers override SetupChangeOfBasisMatrix() which re-sets this for every element and quad point.

Definition at line 258 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver().

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mCheckedOutwardNormals
protected

Whether the boundary elements of the mesh have been checked for whether the ordering if such that the normals are outward-facing

Definition at line 293 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::mCurrentTime
protected

This solver is for static problems, however the body force or surface tractions could be a function of time. The user should call SetCurrentTime() if this is the case.

Definition at line 287 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::SetCurrentTime().

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mFirstStep
protected

Whether this is the first Newton iteration or not

Definition at line 311 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::MIN_NEWTON_ABS_TOL = 1e-10
staticprotected

Minimum absolute tolerance for Newton solve. See documentation for MAX_NEWTON_ABS_TOL.

Definition at line 235 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mIncludeActiveTension
protected

Whether to call AddActiveStressAndStressDerivative() when computing stresses or not.

Subclasses, such as the cardiac mechanics solvers, may implement the above method to add active contributions to the stress. However, sometimes we might want to switch this off. Defaults to true.

Definition at line 331 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::SetIncludeActiveTension().

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::mKspAbsoluteTol
protected

Absolute tolerance for linear systems. Can be set by calling SetKspAbsoluteTolerances(), but default to -1, in which case a relative tolerance is used.

Definition at line 265 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::SetKspAbsoluteTolerance().

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::mLastDampingValue
protected

The last damping value used in the current nonlinear (non-snes) solve. If near 1.0, this indicates that the current guess is near the solution. Initialised to 0.0 at the beginning of each nonlinear solve.

Definition at line 306 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
unsigned AbstractNonlinearElasticitySolver< DIM >::mNumNewtonIterations
protected

Number of Newton iterations taken in last solve.

Definition at line 280 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mPetscDirectSolve
protected
template<unsigned DIM>
Mat& AbstractNonlinearElasticitySolver< DIM >::mrJacobianMatrix
protected

Reference to the matrix 'mLhsSystemMatrix' in the parent class, named as the jacobian, just to make code clearer.

Definition at line 250 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
SolidMechanicsProblemDefinition<DIM>& AbstractNonlinearElasticitySolver< DIM >::mrProblemDefinition
protected

This class contains all the information about the problem (except the material law): body force, surface tractions, fixed nodes, density

Definition at line 244 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver().

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mSetComputeAverageStressPerElement
protected

The user may request that the stress for each element (averaged over quadrature point stresses) are saved during the Solve; this bool states this (defaults to false).

Definition at line 337 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mTakeFullFirstNewtonStep
protected
template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mUseSnesSolver
protected

If SetUsingSnesSolver() is called on the problem definition class, or the command line argument "-mech_use_snes" is given the Petsc SNES solver (nonlinear solver) will be used.

Definition at line 299 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver().

template<unsigned DIM>
bool AbstractNonlinearElasticitySolver< DIM >::mWriteOutputEachNewtonIteration
protected

By default only the initial and final solutions are written. However, we may want to write the solutions after every Newton iteration, in which case the following should be set to true.

Definition at line 272 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::SetWriteOutputEachNewtonIteration().

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::NEWTON_REL_TOL = 1e-4
staticprotected

Relative tolerance for Newton solve. See documentation for MAX_NEWTON_ABS_TOL.

Definition at line 238 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_NODES_PER_BOUNDARY_ELEMENT = DIM*(DIM+1)/2
staticprotected

Number of nodes per boundary element.

Definition at line 217 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_NODES_PER_ELEMENT = (DIM+1)*(DIM+2)/2
staticprotected

Number of nodes per element.

Definition at line 214 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_VERTICES_PER_ELEMENT = DIM+1
staticprotected

Number of vertices per element.

Definition at line 211 of file AbstractNonlinearElasticitySolver.hpp.


The documentation for this class was generated from the following file: