Chaste  Release::2017.1
AbstractContinuumMechanicsSolver< DIM > Class Template Referenceabstract

#include <AbstractContinuumMechanicsSolver.hpp>

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

Public Member Functions

 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 ()
 
virtual std::vector< c_vector< double, DIM > > & rGetSpatialSolution ()=0
 
std::vector< double > & rGetPressures ()
 

Protected Member Functions

void AllocateMatrixMemory ()
 
void ApplyDirichletBoundaryConditions (ApplyDirichletBcsType type, bool symmetricProblem)
 
void AddIdentityBlockForDummyPressureVariables (ApplyDirichletBcsType type)
 
void RemovePressureDummyValuesThroughLinearInterpolation ()
 

Protected Attributes

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
 

Friends

class StressRecoveror< DIM >
 
class VtkNonlinearElasticitySolutionWriter< DIM >
 

Detailed Description

template<unsigned DIM>
class AbstractContinuumMechanicsSolver< DIM >

General base class for continuum mechanics solvers. Deals with memory allocation, writing to file, applying boundary conditions, and adding an identity block (see below)

Note: the PROBLEM DIMENSION depends on DIM and whether a compressible or incompressible problem is being solved. (i) Compressible: mProblemDimension is set to DIM, and the ordering of the unknowns is, for 2D say: [U1 V1 U2 V2 .. Un Vn] (ii) Incompressible: Here we have spatial unknowns (at all nodes of the QuadraticMesh) and pressure unknowns (only at the vertices, as linear bases are). We choose mProblemDim = DIM+1, use (for parallelisation reasons) the ordering [U1 V1 P1 , .. , Un Vn, Pn], introducing dummy variables for pressure unknowns at internal hences, with the equation Pi=0 at these nodes, hence the need for an identity block to be added the corresponding parts of the matrix.

Definition at line 87 of file AbstractContinuumMechanicsSolver.hpp.

Constructor & Destructor Documentation

Member Function Documentation

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::AddIdentityBlockForDummyPressureVariables ( ApplyDirichletBcsType  type)
protected

For incompressible problems, we use the following ordering: [U0 V0 W0 P0 U1 V1 W1 P1 .. Un Vn Wn Pn] where (U V W) is the displacement, and P is the pressure. Therefore P is defined at every node; however for the solve, linear basis functions are used for the pressure, so pressure is solved for at each vertex, not at internal nodes (which are for quadratic basis functions).

The pressure variables for non-vertex nodes are therefore dummy variables. This method enforces the condition P_i=0, where i corresponds to a non-vertex node.

The first input parameter should be one of the following LINEAR_PROBLEM – indicating the overall problem is linear, and in which case the matrix will be altered (1 on diagonal, zeros assumed on rest of row) and the rhs vector will be set to 0.

NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY – indicating the overall problem is nonlinear and here only the residual vector will be altered (sets the residual value of the appropriate rows to P_i-0.0.

NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING – indicating the overall problem is nonlinear, and here, the residual vector will be altered, as will the matrix and RHS vector (see documentation on mResidualVector for why there is a separate residual vector and RHS vector).

Note the row (and columns) for the dummy variables are not explicitly zeroed, as they are assumed to be zero already - in a finite element assembly nothing will have been written to the rows or columns corresponding to dummy variables. Hence this method always maintains symmetry.

Parameters
typesee above

Definition at line 796 of file AbstractContinuumMechanicsSolver.hpp.

References AbstractMesh< ELEMENT_DIM, SPACE_DIM >::GetNode(), AbstractMesh< ELEMENT_DIM, SPACE_DIM >::GetNumNodes(), AbstractContinuumMechanicsSolver< DIM >::mCompressibilityType, AbstractContinuumMechanicsSolver< DIM >::mCurrentSolution, AbstractContinuumMechanicsSolver< DIM >::mLinearSystemRhsVector, AbstractContinuumMechanicsSolver< DIM >::mPreconditionMatrix, AbstractContinuumMechanicsSolver< DIM >::mResidualVector, AbstractContinuumMechanicsSolver< DIM >::mrQuadMesh, AbstractContinuumMechanicsSolver< DIM >::mSystemLhsMatrix, PetscMatTools::SetElement(), and PetscVecTools::SetElement().

Referenced by StokesFlowSolver< DIM >::AssembleSystem().

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::ApplyDirichletBoundaryConditions ( ApplyDirichletBcsType  type,
bool  symmetricProblem 
)
protected

Apply the Dirichlet boundary conditions to the linear system.

The first input parameter should be one of the following LINEAR_PROBLEM – indicating the overall problem is linear, and in which case the BCs will be applied to both the matrix and vector of the linear system

NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY – indicating the overall problem is nonlinear and here only the residual vector will be altered (apply the Dirichlet boundary conditions to the residual vector involves setting appropriate components to the difference between the current value and the correct value).

NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING – indicating the overall problem is nonlinear, and here, the residual vector will be altered, as will the matrix and RHS vector (see documentation on mResidualVector for why there is a separate residual vector and RHS vector).

The second parameter should be true if the overall problem is symmetric, in which case boundary conditions will be applied keeping the matrix symmetric (if the matrix is being altered). See in-code comments for how this is done.

Parameters
typesee above
symmetricProblemsee above

Definition at line 672 of file AbstractContinuumMechanicsSolver.hpp.

References PetscVecTools::AddScaledVector(), PetscTools::Destroy(), PetscMatTools::Finalise(), PetscMatTools::GetMatrixRowDistributed(), AbstractContinuumMechanicsSolver< DIM >::mCurrentSolution, AbstractContinuumMechanicsSolver< DIM >::mDirichletBoundaryConditionsVector, AbstractContinuumMechanicsSolver< DIM >::mLinearSystemRhsVector, AbstractContinuumMechanicsSolver< DIM >::mPreconditionMatrix, AbstractContinuumMechanicsSolver< DIM >::mProblemDimension, AbstractContinuumMechanicsSolver< DIM >::mResidualVector, AbstractContinuumMechanicsSolver< DIM >::mrProblemDefinition, AbstractContinuumMechanicsSolver< DIM >::mSystemLhsMatrix, PetscVecTools::SetElement(), PetscVecTools::Zero(), PetscMatTools::ZeroRowsAndColumnsWithValueOnDiagonal(), and PetscMatTools::ZeroRowsWithValueOnDiagonal().

Referenced by StokesFlowSolver< DIM >::AssembleSystem(), and AbstractNonlinearElasticitySolver< DIM >::FinishAssembleSystem().

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::RemovePressureDummyValuesThroughLinearInterpolation ( )
protected

For incompressible problems, we use the following ordering: [U0 V0 W0 P0 U1 V1 W1 P1 .. Un Vn Wn Pn] where (U V W) is the displacement, and P is the pressure. Therefore P is defined at every node; however for the solve, linear basis functions are used for the pressure, so pressure is solved for at each vertex, not at internal nodes (which are for quadratic basis functions).

The above method, AddIdentityBlockForDummyPressureVariables(), enforces the condition P_i=0, where i corresponds to a non-vertex node.

AFTER the solve, this method can be used to remove the dummy values, but looping over all edges, and linearly interpolating the pressure at the two vertices onto the internal node.

This method assumes each internal node is midway between the two vertices.

Definition at line 593 of file AbstractContinuumMechanicsSolver.hpp.

References AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM >::GetElementIteratorBegin(), AbstractTetrahedralMesh< ELEMENT_DIM, SPACE_DIM >::GetElementIteratorEnd(), AbstractContinuumMechanicsSolver< DIM >::mCurrentSolution, AbstractContinuumMechanicsSolver< DIM >::mProblemDimension, and AbstractContinuumMechanicsSolver< DIM >::mrQuadMesh.

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

template<unsigned DIM>
std::vector<double>& AbstractContinuumMechanicsSolver< DIM >::rGetCurrentSolution ( )
inline
Returns
the current solution vector (advanced use only - for getting the deformed position use rGetDeformedPosition()).

Definition at line 353 of file AbstractContinuumMechanicsSolver.hpp.

References AbstractContinuumMechanicsSolver< DIM >::mCurrentSolution, AbstractContinuumMechanicsSolver< DIM >::rGetPressures(), and AbstractContinuumMechanicsSolver< DIM >::rGetSpatialSolution().

template<unsigned DIM>
std::vector< double > & AbstractContinuumMechanicsSolver< DIM >::rGetPressures ( )
Returns
the pressure, for each NODE in the mesh. If the node is an internal node of the quadratic mesh the pressure is not computed at this node (as linear basis functions are used for the pressure, so pressures unknowns are only present at vertices), so a dummy pressure value of 0 is returned in this vector.

Only valid if mCompressibilityType==INCOMPRESSIBLE

Definition at line 578 of file AbstractContinuumMechanicsSolver.hpp.

References AbstractMesh< ELEMENT_DIM, SPACE_DIM >::GetNumNodes(), AbstractContinuumMechanicsSolver< DIM >::mCurrentSolution, AbstractContinuumMechanicsSolver< DIM >::mPressureSolution, AbstractContinuumMechanicsSolver< DIM >::mProblemDimension, and AbstractContinuumMechanicsSolver< DIM >::mrQuadMesh.

Referenced by AbstractContinuumMechanicsSolver< DIM >::CreateVtkOutput(), AbstractContinuumMechanicsSolver< DIM >::rGetCurrentSolution(), and AbstractContinuumMechanicsSolver< DIM >::WriteCurrentPressureSolution().

template<unsigned DIM>
virtual std::vector<c_vector<double,DIM> >& AbstractContinuumMechanicsSolver< DIM >::rGetSpatialSolution ( )
pure virtual
Returns
the spatial solution. For solids problems this will be the deformed position, for fluids problems this will be the flow.

Implemented in AbstractNonlinearElasticitySolver< DIM >, and StokesFlowSolver< DIM >.

Referenced by AbstractContinuumMechanicsSolver< DIM >::CreateVtkOutput(), AbstractContinuumMechanicsSolver< DIM >::rGetCurrentSolution(), and AbstractContinuumMechanicsSolver< DIM >::WriteCurrentSpatialSolution().

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::SetWriteOutput ( bool  writeOutput = true)

Set whether to write any output.

Parameters
writeOutput(defaults to true)

Definition at line 536 of file AbstractContinuumMechanicsSolver.hpp.

References EXCEPTION, AbstractContinuumMechanicsSolver< DIM >::mOutputDirectory, and AbstractContinuumMechanicsSolver< DIM >::mWriteOutput.

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::WriteCurrentPressureSolution ( int  counterToAppend = -1)

Write the pressure solution. Only valid if mCompressibilityType==INCOMPRESSIBLE. Writes the pressure for ALL nodes on the mesh, including internal nodes (as these are not assumed to be have indices greater than vertex nodes). As linear basis functions are used for pressure, the pressure solution is only computed at the vertices, and pressure dummy variables are used at internal nodes, hence for each internal node, 0 will be written to file.

Parameters
counterToAppendappend a counter to the file name

WriteCurrentPressureSolution() –> file called "pressure.txt" WriteCurrentPressureSolution(3) –> file called "pressure_3.txt"

Definition at line 500 of file AbstractContinuumMechanicsSolver.hpp.

References PetscTools::AmMaster(), PetscTools::Barrier(), AbstractMesh< ELEMENT_DIM, SPACE_DIM >::GetNode(), AbstractContinuumMechanicsSolver< DIM >::mpOutputFileHandler, AbstractContinuumMechanicsSolver< DIM >::mrQuadMesh, AbstractContinuumMechanicsSolver< DIM >::mWriteOutput, OutputFileHandler::OpenOutputFile(), and AbstractContinuumMechanicsSolver< DIM >::rGetPressures().

Referenced by StokesFlowSolver< DIM >::Solve().

template<unsigned DIM>
void AbstractContinuumMechanicsSolver< DIM >::WriteCurrentSpatialSolution ( std::string  fileName,
std::string  fileExtension,
int  counterToAppend = -1 
)

Write the spatial solution (deformed position if solids, flow if fluids) at the nodes

Parameters
fileName(stem)
fileExtensionto append at end.
counterToAppendappend a counter to the file name (defaults to nothing appended).

For example: WriteCurrentSpatialSolution("solution","nodes") –> file called "solution.nodes" WriteCurrentSpatialSolution("solution","nodes",3) –> file called "solution_3.nodes"

Definition at line 458 of file AbstractContinuumMechanicsSolver.hpp.

References PetscTools::AmMaster(), PetscTools::Barrier(), AbstractContinuumMechanicsSolver< DIM >::mpOutputFileHandler, AbstractContinuumMechanicsSolver< DIM >::mWriteOutput, OutputFileHandler::OpenOutputFile(), and AbstractContinuumMechanicsSolver< DIM >::rGetSpatialSolution().

Referenced by StokesFlowSolver< DIM >::Solve(), AbstractNonlinearElasticitySolver< DIM >::Solve(), and AbstractNonlinearElasticitySolver< DIM >::SolveNonSnes().

Member Data Documentation

template<unsigned DIM>
Vec AbstractContinuumMechanicsSolver< DIM >::mDirichletBoundaryConditionsVector
protected
template<unsigned DIM>
std::vector<double> AbstractContinuumMechanicsSolver< DIM >::mPressureSolution
protected

Pressures solution at each vertex of the mesh. Only valid if mCompressibilityType==INCOMPRESSIBLE.

Definition at line 117 of file AbstractContinuumMechanicsSolver.hpp.

Referenced by AbstractContinuumMechanicsSolver< DIM >::rGetPressures().

template<unsigned DIM>
unsigned AbstractContinuumMechanicsSolver< DIM >::mProblemDimension
protected

The problem dimension - the number of unknowns at each node. For compressible problems, where only the deformation/flow is computed at each node, this is equal to DIM For incompressible problems, where pressure is computed as well, this is equal to DIM+1 (there is a pressure variable defined even for internal nodes at which pressure is not computed, this is a dummy pressure variable – done like this for parallelisation reasons. See above).

Definition at line 148 of file AbstractContinuumMechanicsSolver.hpp.

Referenced by AbstractContinuumMechanicsSolver< DIM >::AbstractContinuumMechanicsSolver(), AbstractContinuumMechanicsSolver< DIM >::AllocateMatrixMemory(), AbstractContinuumMechanicsSolver< DIM >::ApplyDirichletBoundaryConditions(), AbstractNonlinearElasticitySolver< DIM >::AssembleOnBoundaryElementForPressureOnDeformedBc(), AbstractNonlinearElasticitySolver< DIM >::GetElementCentroidStrain(), AbstractContinuumMechanicsSolver< DIM >::RemovePressureDummyValuesThroughLinearInterpolation(), AbstractContinuumMechanicsSolver< DIM >::rGetPressures(), and AbstractNonlinearElasticitySolver< DIM >::rGetSpatialSolution().

template<unsigned DIM>
Vec AbstractContinuumMechanicsSolver< DIM >::mResidualVector
protected

Residual vector nonlinear problems.

Since the residual in nonlinear problems is usually also the RHS vector in the linear system, it may seem unncessary to also have the member variable mLinearSystemRhsVector.

However: Newton's method is Ju = f, where J is the Jacobian, u the (negative of the) update and f the residual, but when applying Dirichlet boundary conditions in the compressible case, we alter the rows of the matrix and also alter the columns in order to maintain symmetry. This requires making further changes to the right-hand vector, meaning that it no longer properly represents the residual. Hence, we have to use two vectors.

Overall, this can be represents as

  • compute residual f
  • compute Jacobian J
  • apply BCs to f.
  • alter the linear system from Ju=f to (J*)u=f* which enforces the dirichlet boundary conditions but enforces them symmetrically.

mLinearSystemRhsVector below represents f*.

Definition at line 182 of file AbstractContinuumMechanicsSolver.hpp.

Referenced by AbstractContinuumMechanicsSolver< DIM >::AddIdentityBlockForDummyPressureVariables(), AbstractContinuumMechanicsSolver< DIM >::AllocateMatrixMemory(), AbstractContinuumMechanicsSolver< DIM >::ApplyDirichletBoundaryConditions(), AbstractNonlinearElasticitySolver< DIM >::CalculateResidualNorm(), AbstractNonlinearElasticitySolver< DIM >::ComputeResidual(), AbstractNonlinearElasticitySolver< DIM >::FinishAssembleSystem(), AbstractNonlinearElasticitySolver< DIM >::SolveSnes(), AbstractNonlinearElasticitySolver< DIM >::TakeNewtonStep(), and AbstractContinuumMechanicsSolver< DIM >::~AbstractContinuumMechanicsSolver().

template<unsigned DIM>
ContinuumMechanicsProblemDefinition<DIM>& AbstractContinuumMechanicsSolver< DIM >::mrProblemDefinition
protected
template<unsigned DIM>
AbstractTetrahedralMesh<DIM, DIM>& AbstractContinuumMechanicsSolver< DIM >::mrQuadMesh
protected

The mesh to be solved on. Requires 6 nodes per triangle (or 10 per tetrahedron) as quadratic bases are used.

Definition at line 97 of file AbstractContinuumMechanicsSolver.hpp.

Referenced by AbstractContinuumMechanicsSolver< DIM >::AbstractContinuumMechanicsSolver(), AbstractContinuumMechanicsSolver< DIM >::AddIdentityBlockForDummyPressureVariables(), AbstractNonlinearElasticitySolver< DIM >::AddStressToAverageStressPerElement(), AbstractContinuumMechanicsSolver< DIM >::AllocateMatrixMemory(), AbstractNonlinearElasticitySolver< DIM >::AssembleOnBoundaryElement(), AbstractNonlinearElasticitySolver< DIM >::AssembleOnBoundaryElementForPressureOnDeformedBc(), AbstractNonlinearElasticitySolver< DIM >::CreateCmguiOutput(), AbstractContinuumMechanicsSolver< DIM >::CreateVtkOutput(), AbstractNonlinearElasticitySolver< DIM >::GetAverageStressPerElement(), AbstractNonlinearElasticitySolver< DIM >::GetElementCentroidStrain(), AbstractContinuumMechanicsSolver< DIM >::RemovePressureDummyValuesThroughLinearInterpolation(), AbstractContinuumMechanicsSolver< DIM >::rGetPressures(), StokesFlowSolver< DIM >::rGetSpatialSolution(), AbstractNonlinearElasticitySolver< DIM >::rGetSpatialSolution(), AbstractNonlinearElasticitySolver< DIM >::SetComputeAverageStressPerElementDuringSolve(), AbstractNonlinearElasticitySolver< DIM >::Solve(), StokesFlowSolver< DIM >::StokesFlowSolver(), AbstractNonlinearElasticitySolver< DIM >::WriteCurrentAverageElementStresses(), AbstractContinuumMechanicsSolver< DIM >::WriteCurrentPressureSolution(), and AbstractNonlinearElasticitySolver< DIM >::WriteCurrentStrains().

template<unsigned DIM>
std::vector<c_vector<double,DIM> > AbstractContinuumMechanicsSolver< DIM >::mSpatialSolution
protected

Spatial solution: For solids problems, mSpatialSolution[i](j) = x_j (new position) for node i. For fluids problems, mSpatialSolution[i](j) = u_j (flow) for node i.

Definition at line 114 of file AbstractContinuumMechanicsSolver.hpp.

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

template<unsigned DIM>
bool AbstractContinuumMechanicsSolver< DIM >::mVerbose
protected

If SetVerboseDuringSolve() is called on the problem definition class, or the command line argument "-mech_verbose" or "-mech_very_verbose" is given, than this bool will be set to true and lots of details about each nonlinear solve (including timing breakdowns) will be printed out

Definition at line 160 of file AbstractContinuumMechanicsSolver.hpp.

Referenced by AbstractContinuumMechanicsSolver< DIM >::AbstractContinuumMechanicsSolver(), AbstractNonlinearElasticitySolver< DIM >::PrintLineSearchResult(), StokesFlowSolver< DIM >::Solve(), AbstractNonlinearElasticitySolver< DIM >::SolveNonSnes(), AbstractNonlinearElasticitySolver< DIM >::SolveSnes(), AbstractNonlinearElasticitySolver< DIM >::TakeNewtonStep(), and AbstractNonlinearElasticitySolver< DIM >::UpdateSolutionUsingLineSearch().


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