#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

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

OutputFileHandler *	mpOutputFileHandler

std::vector< c_vector< double, DIM > >	mSpatialSolution

std::vector< double >	mPressureSolution

std::vector< double >	mCurrentSolution

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

rQuadMesh	the quadratic mesh
rProblemDefinition	an object defining in particular the body force and boundary conditions
outputDirectory	output directory
compressibilityType	Should 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 812 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 838 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

rC	The Lagrangian deformation tensor (F^T F)
elementIndex	Index of the current element
currentQuadPointGlobalIndex	The index (assuming an outer loop over elements and an inner loop over quadrature points), of the current quadrature point.
rT	The stress to be added to
rDTdE	the stress derivative to be added to, assuming the final parameter is true
addToDTdE	A boolean flag saying whether the stress derivative is required or not.

Definition at line 420 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

rT	2nd PK stress (matrix is assumed symmetric)
elementIndex	element index

Definition at line 991 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

rBoundaryElement	the boundary element to be integrated on
rAelem	The element's contribution to the LHS matrix is returned. There is no need to zero this matrix before calling.
rBelem	The element's contribution to the RHS vector is returned. There is no need to zero this vector before calling.
assembleResidual	A bool stating whether to assemble the residual vector.
assembleJacobian	A bool stating whether to assemble the Jacobian matrix.
boundaryConditionIndex	index of this boundary (in the vectors in the problem definition object, in which the boundary conditions are stored

Definition at line 1164 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

rBoundaryElement	the boundary element to be integrated on
rAelem	The element's contribution to the LHS matrix is returned. There is no need to zero this matrix before calling.
rBelem	The element's contribution to the RHS vector is returned. There is no need to zero this vector before calling.
assembleResidual	A bool stating whether to assemble the residual vector.
assembleJacobian	A bool stating whether to assemble the Jacobian matrix.
boundaryConditionIndex	index of this boundary (in the vectors in the problem definition object, in which the boundary conditions are stored

Definition at line 1265 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

assembleResidual	A bool stating whether to assemble the residual vector.
assembleLinearSystem	A 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 1678 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

currentGuess	Input, the current guess for the solution
pJacobian	Output, the jacobian matrix at this guess
pPreconditioner	Output, the preconditioner matrix

Definition at line 2244 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

currentGuess	Input, the current guess for the solution
residualVector	Output, the residual vector given this guess.

Definition at line 2228 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

allowException Sometimes 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 1652 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 962 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

assembleResidual	see documentation for AssembleSystem
assembleLinearSystem	see documentation for AssembleSystem

Definition at line 843 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

elementIndex elementIndex

Returns: stress tensor

Definition at line 1024 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

strainType	Which strain to compute, should be one of: DEFORMATION_GRADIENT_F, DEFORMATION_TENSOR_C, or LAGRANGE_STRAIN_E
rElement	The element
rDeformationGradient	Reference to a matrix, which will be filled in by this method.

Definition at line 1066 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 2117 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

counter	current newton iteration number
normResidual	norm of the residual

Definition at line 2027 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

s	s
residNorm	residual norm.

Definition at line 1879 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 893 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 878 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

setComputeAverageStressPerElement whether to compute stresses (defaults to true)

Definition at line 981 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

time	current time

Definition at line 734 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

includeActiveTension whether to include active tension

Definition at line 656 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

kspAbsoluteTolerance the tolerance

Definition at line 685 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

solver KSP solver object (Petsc object)

Todo:: #2057 Make better choices..

Definition at line 1566 of file AbstractNonlinearElasticitySolver.hpp.

References PetscTools::IsSequential(), and PetscTools::SetOption().

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

takeFullFirstStep Whether to take a full first Newton step or not.

Definition at line 704 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

elementIndex	element global index
currentQuadPointGlobalIndex	global index of the quad point

Definition at line 437 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

usePetscDirectSolve Whether to use the Petsc direct solver or not

Definition at line 722 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

writeOutputEachNewtonIteration whether to write each iteration

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

template<unsigned DIM>

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

Solve the problem.

Parameters

tol	tolerance used in Newton solve (defaults to -1.0). Not used in SNES solves.

Definition at line 1522 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

tol	absolute solver used in nonlinear solve

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

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

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

solution The solution of the linear solve in newton's method, ie the update vector u.

Returns: norm of residual

Definition at line 1888 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

rX	X
rY	Y (replicatable vector)
a	a
rZ	Z the returned vector

Definition at line 1688 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

fileName	The file name stem
counterToAppend	(Optional) number to append in the filename.

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

Definition at line 936 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

strainType	Which strain to write, should be one of: DEFORMATION_GRADIENT_F, DEFORMATION_TENSOR_C, or LAGRANGE_STRAIN_E
fileName	The file name stem
counterToAppend	(Optional) number to append in the filename.

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

Definition at line 899 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 217 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 271 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 340 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 226 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 252 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 287 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 281 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 305 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 229 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 325 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 259 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 300 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>

unsigned AbstractNonlinearElasticitySolver< DIM >::mNumNewtonIterations

protected

Number of Newton iterations taken in last solve.

Definition at line 274 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>

bool AbstractNonlinearElasticitySolver< DIM >::mPetscDirectSolve

protected

Get Petsc to do a (almost) direct solve

Definition at line 316 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver(), and AbstractNonlinearElasticitySolver< DIM >::SetUsePetscDirectSolve().

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 244 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 238 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 331 of file AbstractNonlinearElasticitySolver.hpp.

template<unsigned DIM>

bool AbstractNonlinearElasticitySolver< DIM >::mTakeFullFirstNewtonStep

protected

Whether to take a full first Newton step or not - see documentation for SetTakeFullFirstNewtonStep()

Definition at line 311 of file AbstractNonlinearElasticitySolver.hpp.

Referenced by AbstractNonlinearElasticitySolver< DIM >::AbstractNonlinearElasticitySolver(), and AbstractNonlinearElasticitySolver< DIM >::SetTakeFullFirstNewtonStep().

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 293 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 266 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 232 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 211 of file AbstractNonlinearElasticitySolver.hpp.