Chaste Release::3.1
Public Member Functions | Protected Member Functions | Protected Attributes | Static Protected Attributes | Private Member Functions
AbstractNonlinearElasticitySolver< DIM > Class Template Reference

#include <AbstractNonlinearElasticitySolver.hpp>

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

List of all members.

Public Member Functions

void ComputeResidual (Vec currentGuess, Vec residualVector)
void ComputeJacobian (Vec currentGuess, Mat *pJacobian, Mat *pPreconditioner)
 AbstractNonlinearElasticitySolver (QuadraticMesh< 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 SetCurrentTime (double time)
void CreateCmguiOutput ()
void WriteCurrentDeformationGradients (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)

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 GetElementCentroidDeformationGradient (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 Attributes

SolidMechanicsProblemDefinition
< DIM > & 		mrProblemDefinition
Mat & mrJacobianMatrix
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 mIncludeActiveTension
bool mSetComputeAverageStressPerElement
std::vector< c_vector< double,
DIM *(DIM+1)/2 > > 		mAverageStressesPerElement

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 ()

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 115 of file AbstractNonlinearElasticitySolver.hpp.

Constructor & Destructor Documentation

template<unsigned DIM>
AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver ( QuadraticMesh< 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 672 of file AbstractNonlinearElasticitySolver.hpp.

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

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

Destructor.

Definition at line 695 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 
) [inline, protected, virtual]

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 321 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 879 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 1022 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 1124 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 
) [protected, pure 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 CompressibleNonlinearElasticitySolver< DIM >, and IncompressibleNonlinearElasticitySolver< DIM >.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::CalculateResidualNorm ( ) [protected]

Calculate |r|_2 / length(r), where r is the current residual vector.

Definition at line 1530 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 2028 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 2012 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

Definition at line 1504 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 845 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 
) [protected, virtual]

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 701 of file AbstractNonlinearElasticitySolver.hpp.

References PetscMatTools::Finalise(), PetscVecTools::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

Definition at line 913 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::GetElementCentroidDeformationGradient ( Element< DIM, DIM > &  rElement,
c_matrix< double, DIM, DIM > &  rDeformationGradient 
) [protected]

Compute the deformation gradient at the centroid at an element

Parameters:
rElementThe element
rDeformationGradientReference to a matrix, which will be filled in by this method.

Definition at line 956 of file AbstractNonlinearElasticitySolver.hpp.

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

template<unsigned DIM>
unsigned AbstractNonlinearElasticitySolver< DIM >::GetNumNewtonIterations ( )

Get number of Newton iterations taken in last solve.

Definition at line 1924 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::PostNewtonStep ( unsigned  counter,
double  normResidual 
) [protected, virtual]

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 1832 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 1727 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
std::vector< c_vector< double, DIM > > & AbstractNonlinearElasticitySolver< DIM >::rGetDeformedPosition ( )

Get the deformed position. Note: return_value[i](j) = x_j for node i. Just calls rGetSpatialSolution().

Definition at line 752 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
std::vector< c_vector< double, DIM > > & AbstractNonlinearElasticitySolver< DIM >::rGetSpatialSolution ( ) [virtual]

Implemented method, returns the deformed position. Note: return_value[i](j) = x_j for node i.

Implements AbstractContinuumMechanicsSolver< DIM >.

Definition at line 737 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 864 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION, and PetscTools::IsParallel().

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 597 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 556 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 585 of file AbstractNonlinearElasticitySolver.hpp.

References AbstractNonlinearElasticitySolver< DIM >::mKspAbsoluteTol.

template<unsigned DIM>
void AbstractNonlinearElasticitySolver< DIM >::SetKspSolverAndPcType ( KSP  solver) [protected, virtual]

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 1437 of file AbstractNonlinearElasticitySolver.hpp.

References PetscTools::IsSequential().

template<unsigned DIM>
virtual void AbstractNonlinearElasticitySolver< DIM >::SetupChangeOfBasisMatrix ( unsigned  elementIndex,
unsigned  currentQuadPointGlobalIndex 
) [inline, protected, virtual]

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 338 of file AbstractNonlinearElasticitySolver.hpp.

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 574 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

Definition at line 1101 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 1390 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 1838 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 1937 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 1553 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.

Definition at line 1736 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 1538 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 797 of file AbstractNonlinearElasticitySolver.hpp.

References EXCEPTION.

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

Write the deformation gradients for each element (evaluated at the centroids of each element) Each line of the output file corresponds to one element: the DIM*DIM matrix will be written as one line, using the ordering: F00 F01 F02 F10 F11 F12 F20 F21 F22.

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

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

Definition at line 759 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 [static, protected]

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).

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

Definition at line 132 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 186 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 240 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::MAX_NEWTON_ABS_TOL = 1e-7 [static, protected]

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 141 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 167 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 202 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 196 of file AbstractNonlinearElasticitySolver.hpp.

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

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::MIN_NEWTON_ABS_TOL = 1e-10 [static, protected]

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

Definition at line 144 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 225 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 174 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 215 of file AbstractNonlinearElasticitySolver.hpp.

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

Number of Newton iterations taken in last solve.

Definition at line 189 of file AbstractNonlinearElasticitySolver.hpp.

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 159 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

Reimplemented from AbstractContinuumMechanicsSolver< DIM >.

Definition at line 153 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 231 of file AbstractNonlinearElasticitySolver.hpp.

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 208 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 181 of file AbstractNonlinearElasticitySolver.hpp.

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

template<unsigned DIM>
double AbstractNonlinearElasticitySolver< DIM >::NEWTON_REL_TOL = 1e-4 [static, protected]

Relative tolerance for Newton solve. See documentation for MAX_NEWTON_ABS_TOL.

Definition at line 147 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_NODES_PER_BOUNDARY_ELEMENT = DIM*(DIM+1)/2 [static, protected]

Number of nodes per boundary element.

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

Definition at line 126 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_NODES_PER_ELEMENT = (DIM+1)*(DIM+2)/2 [static, protected]

Number of nodes per element.

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

Definition at line 123 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>
const size_t AbstractNonlinearElasticitySolver< DIM >::NUM_VERTICES_PER_ELEMENT = DIM+1 [static, protected]

Number of vertices per element.

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

Definition at line 120 of file AbstractNonlinearElasticitySolver.hpp.

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