36#ifndef ABSTRACTNONLINEARELASTICITYSOLVER_HPP_
37#define ABSTRACTNONLINEARELASTICITYSOLVER_HPP_
41#include "AbstractContinuumMechanicsSolver.hpp"
42#include "LinearSystem.hpp"
44#include "MechanicsEventHandler.hpp"
45#include "ReplicatableVector.hpp"
46#include "FourthOrderTensor.hpp"
47#include "CmguiDeformedSolutionsWriter.hpp"
48#include "AbstractMaterialLaw.hpp"
49#include "QuadraticBasisFunction.hpp"
50#include "SolidMechanicsProblemDefinition.hpp"
52#include "AbstractPerElementWriter.hpp"
61typedef enum StrainType_
63 DEFORMATION_GRADIENT_F = 0,
69#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2)
87PetscErrorCode AbstractNonlinearElasticitySolver_ComputeResidual(SNES snes,
92#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
105 template<
unsigned DIM>
106 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
123template<
unsigned DIM>
124PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
126 Mat* pGlobalJacobian,
127 Mat* pPreconditioner,
128 MatStructure* pMatStructure,
132template <
unsigned DIM>
141template<
unsigned DIM>
171 assert(localIndex == pElement->
GetIndex());
173 c_matrix<double, DIM, DIM> data =
mpSolver->GetAverageStressPerElement(localIndex);
174 for (
unsigned i=0; i<DIM; i++)
176 for (
unsigned j=0; j<DIM; j++)
178 rData[i*DIM+j] = data(i,j);
195template<
unsigned DIM>
400 c_matrix<double,DIM,DIM>& rDeformationGradient);
421 unsigned elementIndex,
422 unsigned currentQuadPointGlobalIndex,
423 c_matrix<double,DIM,DIM>& rT,
460 c_matrix<double, BOUNDARY_STENCIL_SIZE, BOUNDARY_STENCIL_SIZE>& rAelem,
461 c_vector<double, BOUNDARY_STENCIL_SIZE>& rBelem,
462 bool assembleResidual,
463 bool assembleJacobian,
464 unsigned boundaryConditionIndex);
503 c_matrix<double,BOUNDARY_STENCIL_SIZE,BOUNDARY_STENCIL_SIZE>& rAelem,
504 c_vector<double,BOUNDARY_STENCIL_SIZE>& rBelem,
505 bool assembleResidual,
506 bool assembleJacobian,
507 unsigned boundaryConditionIndex);
578 virtual void PostNewtonStep(
unsigned counter,
double normResidual);
632 std::string outputDirectory,
633 CompressibilityType compressibilityType);
645 void Solve(
double tol=-1.0);
687 assert(kspAbsoluteTolerance > 0);
759 void WriteCurrentStrains(StrainType strainType, std::string fileName,
int counterToAppend = -1);
811template<
unsigned DIM>
814 std::string outputDirectory,
815 CompressibilityType compressibilityType)
817 mrProblemDefinition(rProblemDefinition),
818 mrJacobianMatrix(this->mSystemLhsMatrix),
820 mWriteOutputEachNewtonIteration(false),
821 mNumNewtonIterations(0),
823 mCheckedOutwardNormals(false),
824 mLastDampingValue(0.0),
825 mIncludeActiveTension(true),
826 mSetComputeAverageStressPerElement(false)
837template<
unsigned DIM>
842template<
unsigned DIM>
847 if (assembleJacobian)
852 VecCopy(this->mResidualVector, this->mLinearSystemRhsVector);
856 if (assembleJacobian)
858 this->ApplyDirichletBoundaryConditions(NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING, this->mCompressibilityType==COMPRESSIBLE);
860 else if (assembleResidual)
862 this->ApplyDirichletBoundaryConditions(NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY, this->mCompressibilityType==COMPRESSIBLE);
865 if (assembleResidual)
869 if (assembleJacobian)
877template<
unsigned DIM>
880 this->mSpatialSolution.clear();
881 this->mSpatialSolution.resize(this->mrQuadMesh.GetNumNodes(), zero_vector<double>(DIM));
882 for (
unsigned i=0; i<this->mrQuadMesh.GetNumNodes(); i++)
884 for (
unsigned j=0; j<DIM; j++)
886 this->mSpatialSolution[i](j) = this->mrQuadMesh.GetNode(i)->rGetLocation()[j] + this->mCurrentSolution[this->mProblemDimension*i+j];
889 return this->mSpatialSolution;
892template<
unsigned DIM>
895 return rGetSpatialSolution();
898template<
unsigned DIM>
901 if (!this->mWriteOutput)
906 std::stringstream file_name;
907 file_name << fileName;
908 if (counterToAppend >= 0)
910 file_name <<
"_" << counterToAppend;
912 file_name <<
".strain";
914 out_stream p_file = this->mpOutputFileHandler->OpenOutputFile(file_name.str());
916 c_matrix<double,DIM,DIM> strain;
919 iter != this->mrQuadMesh.GetElementIteratorEnd();
922 GetElementCentroidStrain(strainType, *iter, strain);
923 for (
unsigned i=0; i<DIM; i++)
925 for (
unsigned j=0; j<DIM; j++)
927 *p_file << strain(i,j) <<
" ";
935template<
unsigned DIM>
938 if (!this->mWriteOutput)
943 if (!mSetComputeAverageStressPerElement)
945 EXCEPTION(
"Call SetComputeAverageStressPerElementDuringSolve() before solve if calling WriteCurrentAverageElementStresses()");
948 std::stringstream file_name;
949 file_name << fileName;
950 if (counterToAppend >= 0)
952 file_name <<
"_" << counterToAppend;
954 file_name <<
".stress";
955 assert(mAverageStressesPerElement.size()==this->mrQuadMesh.GetNumElements());
958 stress_writer.
WriteData(*(this->mpOutputFileHandler), file_name.str());
961template<
unsigned DIM>
964 if (this->mOutputDirectory ==
"")
966 EXCEPTION(
"No output directory was given so no output was written, cannot convert to cmgui format");
972 WRITE_QUADRATIC_MESH);
974 std::vector<c_vector<double,DIM> >& r_deformed_positions = this->rGetDeformedPosition();
980template<
unsigned DIM>
983 mSetComputeAverageStressPerElement = setComputeAverageStressPerElement;
984 if (setComputeAverageStressPerElement && mAverageStressesPerElement.size()==0)
986 mAverageStressesPerElement.resize(this->mrQuadMesh.GetNumElements(), zero_vector<double>(DIM*(DIM+1)/2));
990template<
unsigned DIM>
993 assert(mSetComputeAverageStressPerElement);
994 assert(elemIndex<this->mrQuadMesh.GetNumElements());
1004 for (
unsigned i=0; i<DIM*(DIM+1)/2; i++)
1015 row = i<=2 ? 0 : (i<=4? 1 : 2);
1016 col = i==0 ? 0 : (i==1 || i==3? 1 : 2);
1019 this->mAverageStressesPerElement[elemIndex](i) += rT(row,col);
1023template<
unsigned DIM>
1026 if (!mSetComputeAverageStressPerElement)
1028 EXCEPTION(
"Call SetComputeAverageStressPerElementDuringSolve() before solve if calling GetAverageStressesPerElement()");
1030 assert(elementIndex<this->mrQuadMesh.GetNumElements());
1032 c_matrix<double,DIM,DIM> stress;
1044 stress(0,0) = mAverageStressesPerElement[elementIndex](0);
1045 stress(1,0) = stress(0,1) = mAverageStressesPerElement[elementIndex](1);
1046 stress(1,1) = mAverageStressesPerElement[elementIndex](2);
1050 stress(0,0) = mAverageStressesPerElement[elementIndex](0);
1051 stress(1,0) = stress(0,1) = mAverageStressesPerElement[elementIndex](1);
1052 stress(2,0) = stress(0,2) = mAverageStressesPerElement[elementIndex](2);
1053 stress(1,1) = mAverageStressesPerElement[elementIndex](3);
1054 stress(2,1) = stress(1,2) = mAverageStressesPerElement[elementIndex](4);
1055 stress(2,2) = mAverageStressesPerElement[elementIndex](5);
1065template<
unsigned DIM>
1068 c_matrix<double,DIM,DIM>& rStrain)
1070 static c_matrix<double,DIM,DIM> jacobian;
1071 static c_matrix<double,DIM,DIM> inverse_jacobian;
1072 double jacobian_determinant;
1074 this->mrQuadMesh.GetInverseJacobianForElement(rElement.
GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
1077 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
1078 for (
unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
1080 for (
unsigned JJ=0; JJ<DIM; JJ++)
1082 element_current_displacements(JJ,II) = this->mCurrentSolution[this->mProblemDimension*rElement.
GetNodeGlobalIndex(II) + JJ];
1087 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
1088 static c_matrix<double,DIM,DIM> grad_u;
1109 grad_u = zero_matrix<double>(DIM,DIM);
1110 for (
unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
1112 for (
unsigned i=0; i<DIM; i++)
1114 for (
unsigned M=0; M<DIM; M++)
1116 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
1121 c_matrix<double,DIM,DIM> deformation_gradient;
1123 for (
unsigned i=0; i<DIM; i++)
1125 for (
unsigned M=0; M<DIM; M++)
1127 deformation_gradient(i,M) = (i==M?1:0) + grad_u(i,M);
1133 case DEFORMATION_GRADIENT_F:
1135 rStrain = deformation_gradient;
1138 case DEFORMATION_TENSOR_C:
1140 rStrain = prod(trans(deformation_gradient),deformation_gradient);
1143 case LAGRANGE_STRAIN_E:
1145 c_matrix<double,DIM,DIM> C = prod(trans(deformation_gradient),deformation_gradient);
1146 for (
unsigned M=0; M<DIM; M++)
1148 for (
unsigned N=0; N<DIM; N++)
1150 rStrain(M,N) = 0.5* ( C(M,N)-(M==N?1:0) );
1163template<
unsigned DIM>
1166 c_matrix<double,BOUNDARY_STENCIL_SIZE,BOUNDARY_STENCIL_SIZE>& rAelem,
1167 c_vector<double,BOUNDARY_STENCIL_SIZE>& rBelem,
1168 bool assembleResidual,
1169 bool assembleJacobian,
1170 unsigned boundaryConditionIndex)
1172 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() == PRESSURE_ON_DEFORMED
1173 || this->mrProblemDefinition.GetTractionBoundaryConditionType() == FUNCTIONAL_PRESSURE_ON_DEFORMED)
1175 AssembleOnBoundaryElementForPressureOnDeformedBc(rBoundaryElement, rAelem, rBelem,
1176 assembleResidual, assembleJacobian, boundaryConditionIndex);
1183 if (assembleJacobian && !assembleResidual)
1189 c_vector<double, DIM> weighted_direction;
1190 double jacobian_determinant;
1191 this->mrQuadMesh.GetWeightedDirectionForBoundaryElement(rBoundaryElement.
GetIndex(), weighted_direction, jacobian_determinant);
1193 c_vector<double,NUM_NODES_PER_BOUNDARY_ELEMENT> phi;
1195 for (
unsigned quad_index=0; quad_index<this->mpBoundaryQuadratureRule->GetNumQuadPoints(); quad_index++)
1197 double wJ = jacobian_determinant * this->mpBoundaryQuadratureRule->GetWeight(quad_index);
1199 const ChastePoint<DIM-1>& quad_point = this->mpBoundaryQuadratureRule->rGetQuadPoint(quad_index);
1205 c_vector<double,DIM> traction = zero_vector<double>(DIM);
1206 switch (this->mrProblemDefinition.GetTractionBoundaryConditionType())
1208 case FUNCTIONAL_TRACTION:
1210 c_vector<double,DIM> X = zero_vector<double>(DIM);
1211 for (
unsigned node_index=0; node_index<NUM_NODES_PER_BOUNDARY_ELEMENT; node_index++)
1213 X += phi(node_index)*this->mrQuadMesh.GetNode( rBoundaryElement.
GetNodeGlobalIndex(node_index) )->rGetLocation();
1215 traction = this->mrProblemDefinition.EvaluateTractionFunction(X, this->mCurrentTime);
1218 case ELEMENTWISE_TRACTION:
1220 traction = this->mrProblemDefinition.rGetElementwiseTractions()[boundaryConditionIndex];
1228 for (
unsigned index=0; index<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index++)
1230 unsigned spatial_dim = index%DIM;
1231 unsigned node_index = (index-spatial_dim)/DIM;
1233 assert(node_index < NUM_NODES_PER_BOUNDARY_ELEMENT);
1235 rBelem(index) -= traction(spatial_dim)
1242template<
unsigned DIM>
1260 return (mLastDampingValue >= 0.5);
1264template<
unsigned DIM>
1267 c_matrix<double,BOUNDARY_STENCIL_SIZE,BOUNDARY_STENCIL_SIZE>& rAelem,
1268 c_vector<double,BOUNDARY_STENCIL_SIZE>& rBelem,
1269 bool assembleResidual,
1270 bool assembleJacobian,
1271 unsigned boundaryConditionIndex)
1273 assert( this->mrProblemDefinition.GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED
1274 || this->mrProblemDefinition.GetTractionBoundaryConditionType()==FUNCTIONAL_PRESSURE_ON_DEFORMED);
1279 c_vector<double, DIM> weighted_direction;
1280 double jacobian_determinant;
1282 this->mrQuadMesh.GetWeightedDirectionForBoundaryElement(rBoundaryElement.
GetIndex(), weighted_direction, jacobian_determinant);
1291 std::set<unsigned> potential_elements = rBoundaryElement.
GetNode(0)->rGetContainingElementIndices();
1292 for (std::set<unsigned>::iterator iter = potential_elements.begin();
1293 iter != potential_elements.end();
1296 p_containing_vol_element = this->mrQuadMesh.GetElement(*iter);
1298 bool this_vol_ele_contains_surf_ele =
true;
1300 for (
unsigned i=1; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
1303 bool found_this_node =
false;
1304 for (
unsigned j=0; j<p_containing_vol_element->
GetNumNodes(); j++)
1307 if (surf_element_node_index == vol_element_node_index)
1309 found_this_node =
true;
1313 if (!found_this_node)
1315 this_vol_ele_contains_surf_ele =
false;
1319 if (this_vol_ele_contains_surf_ele)
1326 std::vector<unsigned> surf_to_vol_map(NUM_NODES_PER_BOUNDARY_ELEMENT);
1327 for (
unsigned i=0; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++)
1330 for (
unsigned j=0; j<NUM_NODES_PER_ELEMENT; j++)
1334 surf_to_vol_map[i] = j;
1343 static c_matrix<double,DIM,DIM> jacobian_vol_element;
1344 static c_matrix<double,DIM,DIM> inverse_jacobian_vol_element;
1345 double jacobian_determinant_vol_element;
1346 this->mrQuadMesh.GetInverseJacobianForElement(p_containing_vol_element->
GetIndex(), jacobian_vol_element, jacobian_determinant_vol_element, inverse_jacobian_vol_element);
1349 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
1350 for (
unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
1352 for (
unsigned JJ=0; JJ<DIM; JJ++)
1354 element_current_displacements(JJ,II) = this->mCurrentSolution[this->mProblemDimension*p_containing_vol_element->
GetNodeGlobalIndex(II) + JJ];
1360 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi_vol_element;
1362 c_vector<double,NUM_NODES_PER_BOUNDARY_ELEMENT> quad_phi_surf_element;
1364 static c_matrix<double, DIM, NUM_NODES_PER_BOUNDARY_ELEMENT> grad_quad_phi_surf_element;
1366 c_matrix<double,DIM,DIM> F;
1367 c_matrix<double,DIM,DIM> invF;
1370 c_matrix<double,1,DIM> normal_as_mat;
1371 for (
unsigned i=0; i<DIM; i++)
1373 normal_as_mat(0,i) = normal(i);
1376 double normal_pressure;
1377 switch (this->mrProblemDefinition.GetTractionBoundaryConditionType())
1379 case PRESSURE_ON_DEFORMED:
1380 normal_pressure = this->mrProblemDefinition.GetNormalPressure();
1382 case FUNCTIONAL_PRESSURE_ON_DEFORMED:
1383 normal_pressure = this->mrProblemDefinition.EvaluateNormalPressureFunction(this->mCurrentTime);
1389 for (
unsigned quad_index=0; quad_index<this->mpBoundaryQuadratureRule->GetNumQuadPoints(); quad_index++)
1391 double wJ = jacobian_determinant * this->mpBoundaryQuadratureRule->GetWeight(quad_index);
1395 const ChastePoint<DIM-1>& quadrature_point = this->mpBoundaryQuadratureRule->rGetQuadPoint(quad_index);
1402 c_vector<double,DIM> X = zero_vector<double>(DIM);
1403 for (
unsigned node_index=0; node_index<NUM_NODES_PER_BOUNDARY_ELEMENT; node_index++)
1405 X += quad_phi_surf_element(node_index)*rBoundaryElement.
GetNode(node_index)->rGetLocation();
1411 c_vector<double,DIM> xi;
1412 for (
unsigned i=0; i<DIM; i++)
1414 xi(i) = weight(i+1);
1420 assert( DIM!=2 || (fabs(weight(0))<1e-6) || (fabs(weight(1))<1e-6) || (fabs(weight(2))<1e-6) );
1424 assert( DIM!=3 || (fabs(weight(0))<1e-6) || (fabs(weight(1))<1e-6) || (fabs(weight(2))<1e-6) || (fabs(weight(3))<1e-6) );
1430 F = identity_matrix<double>(DIM,DIM);
1431 for (
unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
1433 for (
unsigned i=0; i<DIM; i++)
1435 for (
unsigned M=0; M<DIM; M++)
1437 F(i,M) += grad_quad_phi_vol_element(M,node_index)*element_current_displacements(i,node_index);
1445 if (assembleResidual)
1447 c_vector<double,DIM> traction = detF*normal_pressure*prod(trans(invF),normal);
1450 for (
unsigned index=0; index<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index++)
1452 unsigned spatial_dim = index%DIM;
1453 unsigned node_index = (index-spatial_dim)/DIM;
1455 assert(node_index < NUM_NODES_PER_BOUNDARY_ELEMENT);
1457 rBelem(index) -= traction(spatial_dim)
1458 * quad_phi_surf_element(node_index)
1466 if (assembleJacobian && ShouldAssembleMatrixTermForPressureOnDeformedBc())
1468 for (
unsigned II=0; II<NUM_NODES_PER_BOUNDARY_ELEMENT; II++)
1470 for (
unsigned N=0; N<DIM; N++)
1472 grad_quad_phi_surf_element(N,II) = grad_quad_phi_vol_element(N,surf_to_vol_map[II]);
1477 for (
unsigned N=0; N<DIM; N++)
1479 for (
unsigned e=0; e<DIM; e++)
1481 for (
unsigned M=0; M<DIM; M++)
1483 for (
unsigned d=0; d<DIM; d++)
1485 tensor1(N,e,M,d) = invF(N,e)*invF(M,d) - invF(M,e)*invF(N,d);
1493 tensor2.template SetAsContractionOnFirstDimension<DIM>( trans(grad_quad_phi_surf_element), tensor1);
1498 tensor3.template SetAsContractionOnThirdDimension<DIM>( normal_as_mat, tensor2);
1500 for (
unsigned index1=0; index1<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index1++)
1502 unsigned spatial_dim1 = index1%DIM;
1503 unsigned node_index1 = (index1-spatial_dim1)/DIM;
1505 for (
unsigned index2=0; index2<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index2++)
1507 unsigned spatial_dim2 = index2%DIM;
1508 unsigned node_index2 = (index2-spatial_dim2)/DIM;
1510 rAelem(index1,index2) -= normal_pressure
1512 * tensor3(node_index2,spatial_dim2,0,spatial_dim1)
1513 * quad_phi_surf_element(node_index1)
1521template<
unsigned DIM>
1525 mrProblemDefinition.Validate();
1532 if (mrProblemDefinition.GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED && mCheckedOutwardNormals==
false)
1534 this->mrQuadMesh.CheckOutwardNormals();
1535 mCheckedOutwardNormals =
true;
1539 this->WriteCurrentSpatialSolution(
"initial",
"nodes");
1553 if (this->mCompressibilityType==INCOMPRESSIBLE)
1555 this->RemovePressureDummyValuesThroughLinearInterpolation();
1559 this->WriteCurrentSpatialSolution(
"solution",
"nodes");
1565template<
unsigned DIM>
1576 KSPGetPC(solver, &pc);
1578 if (mPetscDirectSolve)
1580 if (this->mCompressibilityType==COMPRESSIBLE)
1582 KSPSetType(solver,KSPPREONLY);
1585 PCSetType(pc, PCLU);
1592 if (this->mCompressibilityType==COMPRESSIBLE)
1594 KSPSetType(solver,KSPCG);
1597 PCSetType(pc, PCICC);
1601 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1)
1609 PCSetType(pc, PCBJACOBI);
1614 unsigned num_restarts = 100;
1615 KSPSetType(solver,KSPGMRES);
1616 KSPGMRESSetRestart(solver,num_restarts);
1618 #ifndef MECH_USE_HYPRE
1619 PCSetType(pc, PCBJACOBI);
1629 PCSetType(pc, PCHYPRE);
1630 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >=2)
1631 KSPSetPCSide(solver, PC_RIGHT);
1633 KSPSetPreconditionerSide(solver, PC_RIGHT);
1651template<
unsigned DIM>
1654 if (!allowException)
1657 AssembleSystem(
true,
false);
1664 AssembleSystem(
true,
false);
1674 return CalculateResidualNorm();
1677template<
unsigned DIM>
1683 VecNorm(this->mResidualVector, NORM_2, &norm);
1684 return norm/this->mNumDofs;
1687template<
unsigned DIM>
1691 std::vector<double>& rZ)
1693 assert(rX.size()==rY.
GetSize());
1694 assert(rY.
GetSize()==rZ.size());
1695 for (
unsigned i=0; i<rX.size(); i++)
1697 rZ[i] = rX[i] + a*rY[i];
1701template<
unsigned DIM>
1713 AssembleSystem(
true,
true);
1726 VecDuplicate(this->mResidualVector,&solution);
1729 KSPCreate(PETSC_COMM_WORLD,&solver);
1731#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
1732 KSPSetOperators(solver, mrJacobianMatrix, this->mPreconditionMatrix);
1734 KSPSetOperators(solver, mrJacobianMatrix, this->mPreconditionMatrix, DIFFERENT_NONZERO_PATTERN );
1738 SetKspSolverAndPcType(solver);
1743 KSPSetFromOptions(solver);
1750 if (mKspAbsoluteTol < 0)
1753 VecDuplicate(this->mResidualVector, &temp);
1755 VecDuplicate(this->mResidualVector, &temp2);
1756 Vec linsys_residual;
1757 VecDuplicate(this->mResidualVector, &linsys_residual);
1759 KSPInitialResidual(solver, solution, temp, temp2, linsys_residual, this->mLinearSystemRhsVector);
1760 double initial_resid_norm;
1761 VecNorm(linsys_residual, NORM_2, &initial_resid_norm);
1767 double ksp_rel_tol = 1e-6;
1768 double absolute_tol = ksp_rel_tol * initial_resid_norm;
1769 if (absolute_tol < 1e-12)
1771 absolute_tol = 1e-12;
1773 KSPSetTolerances(solver, 1e-16, absolute_tol, PETSC_DEFAULT, 1000 );
1777 KSPSetTolerances(solver, 1e-16, mKspAbsoluteTol, PETSC_DEFAULT, 1000 );
1785 KSPSolve(solver,this->mLinearSystemRhsVector,solution);
1811 KSPConvergedReason reason;
1812 KSPGetConvergedReason(solver,&reason);
1813 if (reason == KSP_DIVERGED_ITS || reason == KSP_DIVERGED_BREAKDOWN)
1822 WARN_ONCE_ONLY(
"Linear solve (within a mechanics solve) didn't converge, but this may not stop nonlinear solve converging")
1837 KSPGetIterationNumber(solver, &num_iters);
1842 EXCEPTION(
"KSP Absolute tolerance was too high, linear system wasn't solved - there will be no decrease in Newton residual. Decrease KspAbsoluteTolerance");
1849 std::cout <<
"[" <<
PetscTools::GetMyRank() <<
"]: Num iterations = " << num_iters <<
"\n" << std::flush;
1868 double new_norm_resid = UpdateSolutionUsingLineSearch(solution);
1874 return new_norm_resid;
1877template<
unsigned DIM>
1882 std::cout <<
"\tTesting s = " << s <<
", |f| = " << residNorm <<
"\n" << std::flush;
1886template<
unsigned DIM>
1889 double initial_norm_resid = CalculateResidualNorm();
1892 std::cout <<
"\tInitial |f| [corresponding to s=0] is " << initial_norm_resid <<
"\n" << std::flush;
1897 std::vector<double> old_solution = this->mCurrentSolution;
1899 std::vector<double> damping_values;
1900 for (
unsigned i=10; i>=1; i--)
1902 damping_values.push_back((
double)i/10.0);
1904 damping_values.push_back(0.05);
1905 assert(damping_values.size()==11);
1910 VectorSum(old_solution, update, -damping_values[index], this->mCurrentSolution);
1911 double current_resid_norm = ComputeResidualAndGetNorm(
true);
1912 PrintLineSearchResult(damping_values[index], current_resid_norm);
1917 VectorSum(old_solution, update, -damping_values[index], this->mCurrentSolution);
1918 double next_resid_norm = ComputeResidualAndGetNorm(
true);
1919 PrintLineSearchResult(damping_values[index], next_resid_norm);
1924 while ( (next_resid_norm==DBL_MAX)
1925 || ( (next_resid_norm < current_resid_norm) && (index<damping_values.size()) ) )
1927 current_resid_norm = next_resid_norm;
1930 VectorSum(old_solution, update, -damping_values[index], this->mCurrentSolution);
1931 next_resid_norm = ComputeResidualAndGetNorm(
true);
1932 PrintLineSearchResult(damping_values[index], next_resid_norm);
1937 unsigned best_index;
1939 if (index==damping_values.size() && (next_resid_norm < current_resid_norm))
1949 current_resid_norm = next_resid_norm;
1950 best_index = index-1;
1957 best_index = index-2;
1961 bool full_first_step = mTakeFullFirstNewtonStep && mFirstStep;
1965 if (initial_norm_resid < current_resid_norm && !full_first_step)
1968 EXCEPTION(
"Residual does not appear to decrease in newton direction, quitting");
1973 if (full_first_step)
1977 std::cout <<
"\tTaking full first Newton step...\n";
1984 std::cout <<
"\tChoosing s = " << damping_values[best_index] <<
"\n" << std::flush;
2019 VectorSum(old_solution, update, -damping_values[best_index], this->mCurrentSolution);
2021 mLastDampingValue = damping_values[best_index];
2022 return current_resid_norm;
2025template<
unsigned DIM>
2030template<
unsigned DIM>
2033 mLastDampingValue = 0;
2035 if (mWriteOutputEachNewtonIteration)
2037 this->WriteCurrentSpatialSolution(
"newton_iteration",
"nodes", 0);
2041 double norm_resid = ComputeResidualAndGetNorm(
false);
2044 std::cout <<
"\nNorm of residual is " << norm_resid <<
"\n";
2047 mNumNewtonIterations = 0;
2048 unsigned iteration_number = 1;
2052 tol = NEWTON_REL_TOL*norm_resid;
2055 if (tol > MAX_NEWTON_ABS_TOL)
2057 tol = MAX_NEWTON_ABS_TOL;
2059 if (tol < MIN_NEWTON_ABS_TOL)
2061 tol = MIN_NEWTON_ABS_TOL;
2068 std::cout <<
"Solving with tolerance " << tol <<
"\n";
2071 while (norm_resid > tol)
2075 std::cout <<
"\n-------------------\n"
2076 <<
"Newton iteration " << iteration_number
2077 <<
":\n-------------------\n";
2081 mFirstStep = (iteration_number==1);
2082 norm_resid = TakeNewtonStep();
2086 std::cout <<
"Norm of residual is " << norm_resid <<
"\n";
2089 if (mWriteOutputEachNewtonIteration)
2091 this->WriteCurrentSpatialSolution(
"newton_iteration",
"nodes", iteration_number);
2094 mNumNewtonIterations = iteration_number;
2096 PostNewtonStep(iteration_number,norm_resid);
2099 if (iteration_number==20)
2102 EXCEPTION(
"Not converged after 20 newton iterations, quitting");
2107 if (norm_resid > tol)
2115template<
unsigned DIM>
2118 return mNumNewtonIterations;
2125template<
unsigned DIM>
2130 VecDuplicate(this->mResidualVector, &initial_guess);
2131 double* p_initial_guess;
2132 VecGetArray(initial_guess, &p_initial_guess);
2134 VecGetOwnershipRange(initial_guess, &lo, &hi);
2135 for (
int global_index=lo; global_index<hi; global_index++)
2137 int local_index = global_index - lo;
2138 p_initial_guess[local_index] = this->mCurrentSolution[global_index];
2140 VecRestoreArray(initial_guess, &p_initial_guess);
2143 Vec snes_residual_vec;
2144 VecDuplicate(this->mResidualVector, &snes_residual_vec);
2148 SNESCreate(PETSC_COMM_WORLD, &snes);
2149 SNESSetFunction(snes, snes_residual_vec, &AbstractNonlinearElasticitySolver_ComputeResidual<DIM>,
this);
2150 SNESSetJacobian(snes, mrJacobianMatrix, this->mPreconditionMatrix, &AbstractNonlinearElasticitySolver_ComputeJacobian<DIM>,
this);
2151#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 4)
2152 SNESSetType(snes, SNESNEWTONLS);
2154 SNESSetType(snes, SNESLS);
2156 SNESSetTolerances(snes, 1e-5, 1e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
2158#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 3)
2159 SNESLineSearch linesearch;
2160 SNESGetSNESLineSearch(snes, &linesearch);
2162 SNESLineSearchSetType(linesearch,
"cp");
2163#elif (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 4)
2164 SNESLineSearch linesearch;
2165 SNESGetLineSearch(snes, &linesearch);
2167 SNESLineSearchSetType(linesearch,
"cp");
2170 SNESSetMaxLinearSolveFailures(snes,100);
2173 SNESGetKSP(snes, &ksp);
2175 KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1000 );
2178 SetKspSolverAndPcType(ksp);
2184 SNESSetFromOptions(snes);
2186#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2)
2187 PetscErrorCode err = SNESSolve(snes, initial_guess);
2192 SNESConvergedReason reason;
2193 SNESGetConvergedReason(snes,&reason);
2198 std::stringstream err_stream;
2203 EXCEPTION(
"Nonlinear Solver failed. PETSc error code: "+err_stream.str()+
" .");
2208 std::stringstream reason_stream;
2209 reason_stream << reason;
2213 EXCEPTION(
"Nonlinear Solver did not converge. PETSc reason code: "+reason_stream.str()+
" .");
2218 SNESGetIterationNumber(snes,&num_iters);
2219 mNumNewtonIterations = num_iters;
2226template<
unsigned DIM>
2234 for (
unsigned i=0; i<guess_repl.
GetSize(); i++)
2236 this->mCurrentSolution[i] = guess_repl[i];
2238 AssembleSystem(
true,
false);
2239 VecCopy(this->mResidualVector, residualVector);
2242template<
unsigned DIM>
2252 assert(mrJacobianMatrix==*pJacobian);
2253 assert(this->mPreconditionMatrix==*pPreconditioner);
2257 for (
unsigned i=0; i<guess_repl.
GetSize(); i++)
2259 this->mCurrentSolution[i] = guess_repl[i];
2262 AssembleSystem(
false,
true);
2266template<
unsigned DIM>
2267PetscErrorCode AbstractNonlinearElasticitySolver_ComputeResidual(SNES snes,
2278template<
unsigned DIM>
2279#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
2280 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
2288 p_solver->
ComputeJacobian(currentGuess, &globalJacobian, &preconditioner);
2292 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
2294 Mat* pGlobalJacobian,
2295 Mat* pPreconditioner,
2296 MatStructure* pMatStructure,
2301 p_solver->
ComputeJacobian(currentGuess, pGlobalJacobian, pPreconditioner);
2309template<
unsigned DIM>
2312template<
unsigned DIM>
2315template<
unsigned DIM>
#define EXCEPTION(message)
T Determinant(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
boost::numeric::ublas::c_matrix< T, 1, 1 > Inverse(const boost::numeric::ublas::c_matrix< T, 1, 1 > &rM)
Node< SPACE_DIM > * GetNode(unsigned localIndex) const
unsigned GetNumNodes() const
unsigned GetNodeGlobalIndex(unsigned localIndex) const
unsigned GetIndex() const
c_matrix< double, DIM, DIM > GetAverageStressPerElement(unsigned elementIndex)
virtual void AssembleSystem(bool assembleResidual, bool assembleLinearSystem)=0
static const size_t NUM_NODES_PER_BOUNDARY_ELEMENT
double ComputeResidualAndGetNorm(bool allowException)
void Solve(double tol=-1.0)
static const size_t BOUNDARY_STENCIL_SIZE
bool ShouldAssembleMatrixTermForPressureOnDeformedBc()
void SetCurrentTime(double time)
virtual ~AbstractNonlinearElasticitySolver()
virtual void SetKspSolverAndPcType(KSP solver)
bool mSetComputeAverageStressPerElement
std::vector< c_vector< double, DIM > > & rGetSpatialSolution()
static double NEWTON_REL_TOL
void SetIncludeActiveTension(bool includeActiveTension=true)
bool mIncludeActiveTension
void SolveNonSnes(double tol=-1.0)
unsigned mNumNewtonIterations
virtual void SetupChangeOfBasisMatrix(unsigned elementIndex, unsigned currentQuadPointGlobalIndex)
void SetUsePetscDirectSolve(bool usePetscDirectSolve=true)
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)
static double MIN_NEWTON_ABS_TOL
std::vector< c_vector< double, DIM *(DIM+1)/2 > > mAverageStressesPerElement
unsigned GetNumNewtonIterations()
AbstractNonlinearElasticitySolver(AbstractTetrahedralMesh< DIM, DIM > &rQuadMesh, SolidMechanicsProblemDefinition< DIM > &rProblemDefinition, std::string outputDirectory, CompressibilityType compressibilityType)
void GetElementCentroidStrain(StrainType strainType, Element< DIM, DIM > &rElement, c_matrix< double, DIM, DIM > &rDeformationGradient)
void SetTakeFullFirstNewtonStep(bool takeFullFirstStep=true)
double UpdateSolutionUsingLineSearch(Vec solution)
static const size_t NUM_VERTICES_PER_ELEMENT
bool mCheckedOutwardNormals
virtual void FinishAssembleSystem(bool assembleResidual, bool assembleLinearSystem)
double CalculateResidualNorm()
bool mWriteOutputEachNewtonIteration
void SetComputeAverageStressPerElementDuringSolve(bool setComputeAverageStressPerElement=true)
void ComputeResidual(Vec currentGuess, Vec residualVector)
FourthOrderTensor< DIM, DIM, DIM, DIM > dTdE
std::vector< c_vector< double, DIM > > & rGetDeformedPosition()
static double MAX_NEWTON_ABS_TOL
void VectorSum(std::vector< double > &rX, ReplicatableVector &rY, double a, std::vector< double > &rZ)
void WriteCurrentStrains(StrainType strainType, std::string fileName, int counterToAppend=-1)
void WriteCurrentAverageElementStresses(std::string fileName, int counterToAppend=-1)
static const size_t NUM_NODES_PER_ELEMENT
void PrintLineSearchResult(double s, double residNorm)
void SetKspAbsoluteTolerance(double kspAbsoluteTolerance)
bool mTakeFullFirstNewtonStep
virtual void PostNewtonStep(unsigned counter, double normResidual)
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)
void AddStressToAverageStressPerElement(c_matrix< double, DIM, DIM > &rT, unsigned elementIndex)
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)
void SetWriteOutputEachNewtonIteration(bool writeOutputEachNewtonIteration=true)
void ComputeJacobian(Vec currentGuess, Mat *pJacobian, Mat *pPreconditioner)
c_matrix< double, DIM, DIM > mChangeOfBasisMatrix
SolidMechanicsProblemDefinition< DIM > & mrProblemDefinition
void WriteData(OutputFileHandler &rHandler, const std::string &rFileName)
c_vector< double, SPACE_DIM > CalculateNormal()
c_vector< double, DIM > & rGetLocation()
bool OptionExists(const std::string &rOption)
static CommandLineArguments * Instance()
c_vector< double, SPACE_DIM+1 > CalculateInterpolationWeights(const ChastePoint< SPACE_DIM > &rTestPoint)
static void BeginEvent(unsigned event)
static void EndEvent(unsigned event)
static void ComputeBasisFunctions(const ChastePoint< ELEMENT_DIM > &rPoint, c_vector< double,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)
static void ComputeTransformedBasisFunctionDerivatives(const ChastePoint< ELEMENT_DIM > &rPoint, const c_matrix< double, ELEMENT_DIM, ELEMENT_DIM > &rInverseJacobian, c_matrix< double, ELEMENT_DIM,(ELEMENT_DIM+1) *(ELEMENT_DIM+2)/2 > &rReturnValue)
StressPerElementWriter(AbstractTetrahedralMesh< DIM, DIM > *pMesh, AbstractNonlinearElasticitySolver< DIM > *pSolver)
void Visit(Element< DIM, DIM > *pElement, unsigned localIndex, c_vector< double, DIM *DIM > &rData)
AbstractNonlinearElasticitySolver< DIM > * mpSolver
static void PrintAndReset(std::string message)