Chaste Commit::ca8ccdedf819b6e02855bc0e8e6f50bdecbc5208
AbstractNonlinearElasticitySolver.hpp
1/*
2
3Copyright (c) 2005-2024, University of Oxford.
4All rights reserved.
5
6University of Oxford means the Chancellor, Masters and Scholars of the
7University of Oxford, having an administrative office at Wellington
8Square, Oxford OX1 2JD, UK.
9
10This file is part of Chaste.
11
12Redistribution and use in source and binary forms, with or without
13modification, are permitted provided that the following conditions are met:
14 * Redistributions of source code must retain the above copyright notice,
15 this list of conditions and the following disclaimer.
16 * Redistributions in binary form must reproduce the above copyright notice,
17 this list of conditions and the following disclaimer in the documentation
18 and/or other materials provided with the distribution.
19 * Neither the name of the University of Oxford nor the names of its
20 contributors may be used to endorse or promote products derived from this
21 software without specific prior written permission.
22
23THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33
34*/
35
36#ifndef ABSTRACTNONLINEARELASTICITYSOLVER_HPP_
37#define ABSTRACTNONLINEARELASTICITYSOLVER_HPP_
38
39#include <vector>
40#include <cmath>
41#include "AbstractContinuumMechanicsSolver.hpp"
42#include "LinearSystem.hpp"
43#include "LogFile.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"
51#include "Timer.hpp"
52#include "AbstractPerElementWriter.hpp"
53#include "petscsnes.h"
54
55//#define MECH_USE_HYPRE // uses HYPRE to solve linear systems, requires PETSc to be installed with HYPRE
56
61typedef enum StrainType_
62{
63 DEFORMATION_GRADIENT_F = 0,
64 DEFORMATION_TENSOR_C,
65 LAGRANGE_STRAIN_E
66} StrainType;
67
68// Bizarrely PETSc 2.2 has this, but doesn't put it in the petscksp.h header...
69#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
70extern PetscErrorCode KSPInitialResidual(KSP,Vec,Vec,Vec,Vec,Vec);
71#endif
72
74// Globals functions used by the SNES solver
76
86template<unsigned DIM>
87PetscErrorCode AbstractNonlinearElasticitySolver_ComputeResidual(SNES snes,
88 Vec currentGuess,
89 Vec residualVector,
90 void* pContext);
91
92#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
105 template<unsigned DIM>
106 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
107 Vec currentGuess,
108 Mat globalJacobian,
109 Mat preconditioner,
110 void* pContext);
111#else
123template<unsigned DIM>
124PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
125 Vec currentGuess,
126 Mat* pGlobalJacobian,
127 Mat* pPreconditioner,
128 MatStructure* pMatStructure,
129 void* pContext);
130#endif
131
132template <unsigned DIM>
133class AbstractNonlinearElasticitySolver; //Forward declaration
134
141template<unsigned DIM>
142class StressPerElementWriter : public AbstractPerElementWriter<DIM, DIM, DIM*DIM>
143{
144private:
146public:
147
154 : AbstractPerElementWriter<DIM, DIM, DIM*DIM>(pMesh),
155 mpSolver(pSolver)
156 {
157 }
158
166 void Visit(Element<DIM, DIM>* pElement, unsigned localIndex, c_vector<double, DIM*DIM>& rData)
167 {
171 assert(localIndex == pElement->GetIndex()); //Will fail when we move to DistributedQuadraticMesh
172 //Flatten the matrix
173 c_matrix<double, DIM, DIM> data = mpSolver->GetAverageStressPerElement(localIndex);
174 for (unsigned i=0; i<DIM; i++)
175 {
176 for (unsigned j=0; j<DIM; j++)
177 {
178 rData[i*DIM+j] = data(i,j);
179 }
180 }
181 }
182};
183
195template<unsigned DIM>
197{
198 friend class StressRecoveror<DIM>;
199 friend class VtkNonlinearElasticitySolutionWriter<DIM>;
200
201
202protected:
203
205 static const size_t NUM_VERTICES_PER_ELEMENT = DIM+1;
206
208 static const size_t NUM_NODES_PER_ELEMENT = (DIM+1)*(DIM+2)/2; // assuming quadratic
209
211 static const size_t NUM_NODES_PER_BOUNDARY_ELEMENT = DIM*(DIM+1)/2; // assuming quadratic
212
218
226 static double MAX_NEWTON_ABS_TOL;
227
229 static double MIN_NEWTON_ABS_TOL;
230
232 static double NEWTON_REL_TOL;
233
239
245
252 c_matrix<double,DIM,DIM> mChangeOfBasisMatrix;
253
260
267
272
275
282
283
288
294
301
306
312
317
326
332
340 std::vector<c_vector<double,DIM*(DIM+1)/2> > mAverageStressesPerElement;
341
349 void AddStressToAverageStressPerElement(c_matrix<double,DIM,DIM>& rT, unsigned elementIndex);
350
359 virtual void SetKspSolverAndPcType(KSP solver);
360
361
374 virtual void AssembleSystem(bool assembleResidual, bool assembleLinearSystem)=0;
375
383 virtual void FinishAssembleSystem(bool assembleResidual, bool assembleLinearSystem);
384
386 //
387 // Element level methods
388 //
390
398 void GetElementCentroidStrain(StrainType strainType,
399 Element<DIM,DIM>& rElement,
400 c_matrix<double,DIM,DIM>& rDeformationGradient);
401
420 virtual void AddActiveStressAndStressDerivative(c_matrix<double,DIM,DIM>& rC,
421 unsigned elementIndex,
422 unsigned currentQuadPointGlobalIndex,
423 c_matrix<double,DIM,DIM>& rT,
425 bool addToDTdE)
426 {
427 // does nothing - subclass needs to overload this if there are active stresses
428 }
429
437 virtual void SetupChangeOfBasisMatrix(unsigned elementIndex, unsigned currentQuadPointGlobalIndex)
438 {
439 }
440
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);
465
466
482
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);
508
510 //
511 // These methods form the non-SNES nonlinear solver
512 //
514
526 double ComputeResidualAndGetNorm(bool allowException);
527
529 double CalculateResidualNorm();
530
540 void VectorSum(std::vector<double>& rX, ReplicatableVector& rY, double a, std::vector<double>& rZ);
541
549 void PrintLineSearchResult(double s, double residNorm);
550
558 double TakeNewtonStep();
559
569 double UpdateSolutionUsingLineSearch(Vec solution);
570
578 virtual void PostNewtonStep(unsigned counter, double normResidual);
579
585 void SolveNonSnes(double tol=-1.0);
586
587
589 //
590 // These methods form the SNES nonlinear solver
591 //
593public: // need to be public as are called by global functions
600 void ComputeResidual(Vec currentGuess, Vec residualVector);
601
609 void ComputeJacobian(Vec currentGuess, Mat* pJacobian, Mat* pPreconditioner);
610
611private:
616 void SolveSnes();
617
618public:
619
631 SolidMechanicsProblemDefinition<DIM>& rProblemDefinition,
632 std::string outputDirectory,
633 CompressibilityType compressibilityType);
634
639
645 void Solve(double tol=-1.0);
646
656 void SetIncludeActiveTension(bool includeActiveTension = true)
657 {
658 mIncludeActiveTension = includeActiveTension;
659 }
660
664 unsigned GetNumNewtonIterations();
665
666
674 void SetWriteOutputEachNewtonIteration(bool writeOutputEachNewtonIteration=true)
675 {
676 mWriteOutputEachNewtonIteration = writeOutputEachNewtonIteration;
677 }
678
685 void SetKspAbsoluteTolerance(double kspAbsoluteTolerance)
686 {
687 assert(kspAbsoluteTolerance > 0);
688 mKspAbsoluteTol = kspAbsoluteTolerance;
689 }
690
704 void SetTakeFullFirstNewtonStep(bool takeFullFirstStep = true)
705 {
706 mTakeFullFirstNewtonStep = takeFullFirstStep;
707 }
708
722 void SetUsePetscDirectSolve(bool usePetscDirectSolve = true)
723 {
724 mPetscDirectSolve = usePetscDirectSolve;
725 }
726
727
734 void SetCurrentTime(double time)
735 {
736 mCurrentTime = time;
737 }
738
743 void CreateCmguiOutput();
744
745
759 void WriteCurrentStrains(StrainType strainType, std::string fileName, int counterToAppend = -1);
760
767 void SetComputeAverageStressPerElementDuringSolve(bool setComputeAverageStressPerElement = true);
768
782 void WriteCurrentAverageElementStresses(std::string fileName, int counterToAppend = -1);
783
788 std::vector<c_vector<double,DIM> >& rGetSpatialSolution();
789
794 std::vector<c_vector<double,DIM> >& rGetDeformedPosition();
795
804 c_matrix<double,DIM,DIM> GetAverageStressPerElement(unsigned elementIndex);
805};
806
808// Implementation: first, the non-nonlinear-solve methods
810
811template<unsigned DIM>
813 SolidMechanicsProblemDefinition<DIM>& rProblemDefinition,
814 std::string outputDirectory,
815 CompressibilityType compressibilityType)
816 : AbstractContinuumMechanicsSolver<DIM>(rQuadMesh, rProblemDefinition, outputDirectory, compressibilityType),
817 mrProblemDefinition(rProblemDefinition),
818 mrJacobianMatrix(this->mSystemLhsMatrix),
819 mKspAbsoluteTol(-1),
820 mWriteOutputEachNewtonIteration(false),
821 mNumNewtonIterations(0),
822 mCurrentTime(0.0),
823 mCheckedOutwardNormals(false),
824 mLastDampingValue(0.0),
825 mIncludeActiveTension(true),
826 mSetComputeAverageStressPerElement(false)
827{
828 mUseSnesSolver = (mrProblemDefinition.GetSolveUsingSnes() ||
829 CommandLineArguments::Instance()->OptionExists("-mech_use_snes") );
830
831 mChangeOfBasisMatrix = identity_matrix<double>(DIM,DIM);
832
833 mTakeFullFirstNewtonStep = CommandLineArguments::Instance()->OptionExists("-mech_full_first_newton_step");
834 mPetscDirectSolve = CommandLineArguments::Instance()->OptionExists("-mech_petsc_direct_solve");
835}
836
837template<unsigned DIM>
841
842template<unsigned DIM>
843void AbstractNonlinearElasticitySolver<DIM>::FinishAssembleSystem(bool assembleResidual, bool assembleJacobian)
844{
845 PetscVecTools::Finalise(this->mResidualVector);
846
847 if (assembleJacobian)
848 {
849 PetscMatTools::SwitchWriteMode(mrJacobianMatrix);
850 PetscMatTools::SwitchWriteMode(this->mPreconditionMatrix);
851
852 VecCopy(this->mResidualVector, this->mLinearSystemRhsVector);
853 }
854
855 // Apply Dirichlet boundary conditions
856 if (assembleJacobian)
857 {
858 this->ApplyDirichletBoundaryConditions(NONLINEAR_PROBLEM_APPLY_TO_EVERYTHING, this->mCompressibilityType==COMPRESSIBLE);
859 }
860 else if (assembleResidual)
861 {
862 this->ApplyDirichletBoundaryConditions(NONLINEAR_PROBLEM_APPLY_TO_RESIDUAL_ONLY, this->mCompressibilityType==COMPRESSIBLE);
863 }
864
865 if (assembleResidual)
866 {
867 PetscVecTools::Finalise(this->mResidualVector);
868 }
869 if (assembleJacobian)
870 {
871 PetscMatTools::Finalise(mrJacobianMatrix);
872 PetscMatTools::Finalise(this->mPreconditionMatrix);
873 PetscVecTools::Finalise(this->mLinearSystemRhsVector);
874 }
875}
876
877template<unsigned DIM>
879{
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++)
883 {
884 for (unsigned j=0; j<DIM; j++)
885 {
886 this->mSpatialSolution[i](j) = this->mrQuadMesh.GetNode(i)->rGetLocation()[j] + this->mCurrentSolution[this->mProblemDimension*i+j];
887 }
888 }
889 return this->mSpatialSolution;
890}
891
892template<unsigned DIM>
894{
895 return rGetSpatialSolution();
896}
897
898template<unsigned DIM>
899void AbstractNonlinearElasticitySolver<DIM>::WriteCurrentStrains(StrainType strainType, std::string fileName, int counterToAppend)
900{
901 if (!this->mWriteOutput)
902 {
903 return;
904 }
905
906 std::stringstream file_name;
907 file_name << fileName;
908 if (counterToAppend >= 0)
909 {
910 file_name << "_" << counterToAppend;
911 }
912 file_name << ".strain";
913
914 out_stream p_file = this->mpOutputFileHandler->OpenOutputFile(file_name.str());
915
916 c_matrix<double,DIM,DIM> strain;
917
918 for (typename AbstractTetrahedralMesh<DIM,DIM>::ElementIterator iter = this->mrQuadMesh.GetElementIteratorBegin();
919 iter != this->mrQuadMesh.GetElementIteratorEnd();
920 ++iter)
921 {
922 GetElementCentroidStrain(strainType, *iter, strain);
923 for (unsigned i=0; i<DIM; i++)
924 {
925 for (unsigned j=0; j<DIM; j++)
926 {
927 *p_file << strain(i,j) << " ";
928 }
929 }
930 *p_file << "\n";
931 }
932 p_file->close();
933}
934
935template<unsigned DIM>
937{
938 if (!this->mWriteOutput)
939 {
940 return;
941 }
942
943 if (!mSetComputeAverageStressPerElement)
944 {
945 EXCEPTION("Call SetComputeAverageStressPerElementDuringSolve() before solve if calling WriteCurrentAverageElementStresses()");
946 }
947
948 std::stringstream file_name;
949 file_name << fileName;
950 if (counterToAppend >= 0)
951 {
952 file_name << "_" << counterToAppend;
953 }
954 file_name << ".stress";
955 assert(mAverageStressesPerElement.size()==this->mrQuadMesh.GetNumElements());
956
957 StressPerElementWriter<DIM> stress_writer(&(this->mrQuadMesh),this);
958 stress_writer.WriteData(*(this->mpOutputFileHandler), file_name.str());
959}
960
961template<unsigned DIM>
963{
964 if (this->mOutputDirectory == "")
965 {
966 EXCEPTION("No output directory was given so no output was written, cannot convert to cmgui format");
967 }
968
969 CmguiDeformedSolutionsWriter<DIM> writer(this->mOutputDirectory + "/cmgui",
970 "solution",
971 this->mrQuadMesh,
972 WRITE_QUADRATIC_MESH);
973
974 std::vector<c_vector<double,DIM> >& r_deformed_positions = this->rGetDeformedPosition();
975 writer.WriteInitialMesh(); // this writes solution_0.exnode and .exelem
976 writer.WriteDeformationPositions(r_deformed_positions, 1); // this writes the final solution as solution_1.exnode
977 writer.WriteCmguiScript(); // writes LoadSolutions.com
978}
979
980template<unsigned DIM>
982{
983 mSetComputeAverageStressPerElement = setComputeAverageStressPerElement;
984 if (setComputeAverageStressPerElement && mAverageStressesPerElement.size()==0)
985 {
986 mAverageStressesPerElement.resize(this->mrQuadMesh.GetNumElements(), zero_vector<double>(DIM*(DIM+1)/2));
987 }
988}
989
990template<unsigned DIM>
991void AbstractNonlinearElasticitySolver<DIM>::AddStressToAverageStressPerElement(c_matrix<double,DIM,DIM>& rT, unsigned elemIndex)
992{
993 assert(mSetComputeAverageStressPerElement);
994 assert(elemIndex<this->mrQuadMesh.GetNumElements());
995
996 // In 2d the matrix is
997 // [T00 T01]
998 // [T10 T11]
999 // where T01 = T10. We store this as a vector
1000 // [T00 T01 T11]
1001 //
1002 // Similarly, for 3d we store
1003 // [T00 T01 T02 T11 T12 T22]
1004 for (unsigned i=0; i<DIM*(DIM+1)/2; i++)
1005 {
1006 unsigned row;
1007 unsigned col;
1008 if (DIM == 2)
1009 {
1010 row = i<=1 ? 0 : 1;
1011 col = i==0 ? 0 : 1;
1012 }
1013 else // DIM == 3
1014 {
1015 row = i<=2 ? 0 : (i<=4? 1 : 2);
1016 col = i==0 ? 0 : (i==1 || i==3? 1 : 2);
1017 }
1018
1019 this->mAverageStressesPerElement[elemIndex](i) += rT(row,col);
1020 }
1021}
1022
1023template<unsigned DIM>
1024c_matrix<double,DIM,DIM> AbstractNonlinearElasticitySolver<DIM>::GetAverageStressPerElement(unsigned elementIndex)
1025{
1026 if (!mSetComputeAverageStressPerElement)
1027 {
1028 EXCEPTION("Call SetComputeAverageStressPerElementDuringSolve() before solve if calling GetAverageStressesPerElement()");
1029 }
1030 assert(elementIndex<this->mrQuadMesh.GetNumElements());
1031
1032 c_matrix<double,DIM,DIM> stress;
1033
1034 // In 2d the matrix is
1035 // [T00 T01]
1036 // [T10 T11]
1037 // where T01 = T10, and was stored as
1038 // [T00 T01 T11]
1039 //
1040 // Similarly, for 3d the matrix was stored as
1041 // [T00 T01 T02 T11 T12 T22]
1042 if (DIM == 2)
1043 {
1044 stress(0,0) = mAverageStressesPerElement[elementIndex](0);
1045 stress(1,0) = stress(0,1) = mAverageStressesPerElement[elementIndex](1);
1046 stress(1,1) = mAverageStressesPerElement[elementIndex](2);
1047 }
1048 else
1049 {
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);
1056 }
1057
1058 return stress;
1059}
1060
1062// Methods at the 'element level'.
1064
1065template<unsigned DIM>
1067 Element<DIM,DIM>& rElement,
1068 c_matrix<double,DIM,DIM>& rStrain)
1069{
1070 static c_matrix<double,DIM,DIM> jacobian;
1071 static c_matrix<double,DIM,DIM> inverse_jacobian;
1072 double jacobian_determinant;
1073
1074 this->mrQuadMesh.GetInverseJacobianForElement(rElement.GetIndex(), jacobian, jacobian_determinant, inverse_jacobian);
1075
1076 // Get the current displacement at the nodes
1077 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
1078 for (unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
1079 {
1080 for (unsigned JJ=0; JJ<DIM; JJ++)
1081 {
1082 element_current_displacements(JJ,II) = this->mCurrentSolution[this->mProblemDimension*rElement.GetNodeGlobalIndex(II) + JJ];
1083 }
1084 }
1085
1086 // Allocate memory for the basis functions values and derivative values
1087 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi;
1088 static c_matrix<double,DIM,DIM> grad_u; // grad_u = (du_i/dX_M)
1089
1090 // we need the point in the canonical element which corresponds to the centroid of the
1091 // version of the element in physical space. This point can be shown to be (1/3,1/3).
1092 ChastePoint<DIM> quadrature_point;
1093 if (DIM == 2)
1094 {
1095 quadrature_point.rGetLocation()(0) = 1.0/3.0;
1096 quadrature_point.rGetLocation()(1) = 1.0/3.0;
1097 }
1098 else
1099 {
1100 assert(DIM==3);
1101 quadrature_point.rGetLocation()(0) = 1.0/4.0;
1102 quadrature_point.rGetLocation()(1) = 1.0/4.0;
1103 quadrature_point.rGetLocation()(2) = 1.0/4.0;
1104 }
1105
1106 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(quadrature_point, inverse_jacobian, grad_quad_phi);
1107
1108 // Interpolate grad_u
1109 grad_u = zero_matrix<double>(DIM,DIM);
1110 for (unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
1111 {
1112 for (unsigned i=0; i<DIM; i++)
1113 {
1114 for (unsigned M=0; M<DIM; M++)
1115 {
1116 grad_u(i,M) += grad_quad_phi(M,node_index)*element_current_displacements(i,node_index);
1117 }
1118 }
1119 }
1120
1121 c_matrix<double,DIM,DIM> deformation_gradient;
1122
1123 for (unsigned i=0; i<DIM; i++)
1124 {
1125 for (unsigned M=0; M<DIM; M++)
1126 {
1127 deformation_gradient(i,M) = (i==M?1:0) + grad_u(i,M);
1128 }
1129 }
1130
1131 switch(strainType)
1132 {
1133 case DEFORMATION_GRADIENT_F:
1134 {
1135 rStrain = deformation_gradient;
1136 break;
1137 }
1138 case DEFORMATION_TENSOR_C:
1139 {
1140 rStrain = prod(trans(deformation_gradient),deformation_gradient);
1141 break;
1142 }
1143 case LAGRANGE_STRAIN_E:
1144 {
1145 c_matrix<double,DIM,DIM> C = prod(trans(deformation_gradient),deformation_gradient);
1146 for (unsigned M=0; M<DIM; M++)
1147 {
1148 for (unsigned N=0; N<DIM; N++)
1149 {
1150 rStrain(M,N) = 0.5* ( C(M,N)-(M==N?1:0) );
1151 }
1152 }
1153 break;
1154 }
1155 default:
1156 {
1158 break;
1159 }
1160 }
1161}
1162
1163template<unsigned DIM>
1165 BoundaryElement<DIM-1,DIM>& rBoundaryElement,
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)
1171{
1172 if (this->mrProblemDefinition.GetTractionBoundaryConditionType() == PRESSURE_ON_DEFORMED
1173 || this->mrProblemDefinition.GetTractionBoundaryConditionType() == FUNCTIONAL_PRESSURE_ON_DEFORMED)
1174 {
1175 AssembleOnBoundaryElementForPressureOnDeformedBc(rBoundaryElement, rAelem, rBelem,
1176 assembleResidual, assembleJacobian, boundaryConditionIndex);
1177 return;
1178 }
1179
1180 rAelem.clear();
1181 rBelem.clear();
1182
1183 if (assembleJacobian && !assembleResidual)
1184 {
1185 // Nothing to do
1186 return;
1187 }
1188
1189 c_vector<double, DIM> weighted_direction;
1190 double jacobian_determinant;
1191 this->mrQuadMesh.GetWeightedDirectionForBoundaryElement(rBoundaryElement.GetIndex(), weighted_direction, jacobian_determinant);
1192
1193 c_vector<double,NUM_NODES_PER_BOUNDARY_ELEMENT> phi;
1194
1195 for (unsigned quad_index=0; quad_index<this->mpBoundaryQuadratureRule->GetNumQuadPoints(); quad_index++)
1196 {
1197 double wJ = jacobian_determinant * this->mpBoundaryQuadratureRule->GetWeight(quad_index);
1198
1199 const ChastePoint<DIM-1>& quad_point = this->mpBoundaryQuadratureRule->rGetQuadPoint(quad_index);
1200
1202
1203 // Get the required traction, interpolating X (slightly inefficiently,
1204 // as interpolating using quad bases) if necessary
1205 c_vector<double,DIM> traction = zero_vector<double>(DIM);
1206 switch (this->mrProblemDefinition.GetTractionBoundaryConditionType())
1207 {
1208 case FUNCTIONAL_TRACTION:
1209 {
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++)
1212 {
1213 X += phi(node_index)*this->mrQuadMesh.GetNode( rBoundaryElement.GetNodeGlobalIndex(node_index) )->rGetLocation();
1214 }
1215 traction = this->mrProblemDefinition.EvaluateTractionFunction(X, this->mCurrentTime);
1216 break;
1217 }
1218 case ELEMENTWISE_TRACTION:
1219 {
1220 traction = this->mrProblemDefinition.rGetElementwiseTractions()[boundaryConditionIndex];
1221 break;
1222 }
1223 default:
1225 }
1226
1227
1228 for (unsigned index=0; index<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index++)
1229 {
1230 unsigned spatial_dim = index%DIM;
1231 unsigned node_index = (index-spatial_dim)/DIM;
1232
1233 assert(node_index < NUM_NODES_PER_BOUNDARY_ELEMENT);
1234
1235 rBelem(index) -= traction(spatial_dim)
1236 * phi(node_index)
1237 * wJ;
1238 }
1239 }
1240}
1241
1242template<unsigned DIM>
1244{
1245 if (mUseSnesSolver)
1246 {
1247 // although not using this in the first few steps might be useful when the deformation
1248 // is large, the snes solver is more robust, so we have this on all the time. (Also because
1249 // for cardiac problems and in timesteps after the initial large deformation, we want this on
1250 // in the first step
1251 return true;
1252
1253 // could do something like this, if make the snes a member variable
1254 //PetscInt iteration_number;
1255 //SNESGetIterationNumber(mSnes,&iteration_number);
1256 //return (iteration_number >= 3);
1257 }
1258 else
1259 {
1260 return (mLastDampingValue >= 0.5);
1261 }
1262}
1263
1264template<unsigned DIM>
1266 BoundaryElement<DIM-1,DIM>& rBoundaryElement,
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)
1272{
1273 assert( this->mrProblemDefinition.GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED
1274 || this->mrProblemDefinition.GetTractionBoundaryConditionType()==FUNCTIONAL_PRESSURE_ON_DEFORMED);
1275
1276 rAelem.clear();
1277 rBelem.clear();
1278
1279 c_vector<double, DIM> weighted_direction;
1280 double jacobian_determinant;
1281 // note: jacobian determinant may be over-written below
1282 this->mrQuadMesh.GetWeightedDirectionForBoundaryElement(rBoundaryElement.GetIndex(), weighted_direction, jacobian_determinant);
1283
1285 // Find the volume element of the mesh which
1286 // contains this boundary element
1288
1289 Element<DIM,DIM>* p_containing_vol_element = nullptr;
1290
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();
1294 iter++)
1295 {
1296 p_containing_vol_element = this->mrQuadMesh.GetElement(*iter);
1297
1298 bool this_vol_ele_contains_surf_ele = true;
1299 // loop over the nodes of boundary element and see if they are in the volume element
1300 for (unsigned i=1; i<NUM_NODES_PER_BOUNDARY_ELEMENT; i++) // don't need to start at 0, given looping over contain elems of node 0
1301 {
1302 unsigned surf_element_node_index = rBoundaryElement.GetNodeGlobalIndex(i);
1303 bool found_this_node = false;
1304 for (unsigned j=0; j<p_containing_vol_element->GetNumNodes(); j++)
1305 {
1306 unsigned vol_element_node_index = p_containing_vol_element->GetNodeGlobalIndex(j);
1307 if (surf_element_node_index == vol_element_node_index)
1308 {
1309 found_this_node = true;
1310 break;
1311 }
1312 }
1313 if (!found_this_node)
1314 {
1315 this_vol_ele_contains_surf_ele = false;
1316 break;
1317 }
1318 }
1319 if (this_vol_ele_contains_surf_ele)
1320 {
1321 break;
1322 }
1323 }
1324
1325 // Find the local node index in the volume element for each node in the boundary element
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++)
1328 {
1329 unsigned index = rBoundaryElement.GetNodeGlobalIndex(i);
1330 for (unsigned j=0; j<NUM_NODES_PER_ELEMENT; j++)
1331 {
1332 if (p_containing_vol_element->GetNodeGlobalIndex(j)==index)
1333 {
1334 surf_to_vol_map[i] = j;
1335 break;
1336 }
1337 }
1338 }
1339
1340
1341 // We require the volume element to compute F, which requires grad_phi on the volume element. For this we will
1342 // need the inverse jacobian for the volume element
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);
1347
1348 // Get the current displacements at each node of the volume element, to be used in computing F
1349 static c_matrix<double,DIM,NUM_NODES_PER_ELEMENT> element_current_displacements;
1350 for (unsigned II=0; II<NUM_NODES_PER_ELEMENT; II++)
1351 {
1352 for (unsigned JJ=0; JJ<DIM; JJ++)
1353 {
1354 element_current_displacements(JJ,II) = this->mCurrentSolution[this->mProblemDimension*p_containing_vol_element->GetNodeGlobalIndex(II) + JJ];
1355 }
1356 }
1357
1358
1359 // We will need both {grad phi_i} for the quadratic bases of the volume element, for computing F..
1360 static c_matrix<double, DIM, NUM_NODES_PER_ELEMENT> grad_quad_phi_vol_element;
1361 // ..the phi_i for each of the quadratic bases of the surface element, for the standard FE assembly part.
1362 c_vector<double,NUM_NODES_PER_BOUNDARY_ELEMENT> quad_phi_surf_element;
1363 // We need this too, which is obtained by taking a subset of grad_quad_phi_vol_element
1364 static c_matrix<double, DIM, NUM_NODES_PER_BOUNDARY_ELEMENT> grad_quad_phi_surf_element;
1365
1366 c_matrix<double,DIM,DIM> F;
1367 c_matrix<double,DIM,DIM> invF;
1368
1369 c_vector<double,DIM> normal = rBoundaryElement.CalculateNormal();
1370 c_matrix<double,1,DIM> normal_as_mat;
1371 for (unsigned i=0; i<DIM; i++)
1372 {
1373 normal_as_mat(0,i) = normal(i);
1374 }
1375
1376 double normal_pressure;
1377 switch (this->mrProblemDefinition.GetTractionBoundaryConditionType())
1378 {
1379 case PRESSURE_ON_DEFORMED:
1380 normal_pressure = this->mrProblemDefinition.GetNormalPressure();
1381 break;
1382 case FUNCTIONAL_PRESSURE_ON_DEFORMED:
1383 normal_pressure = this->mrProblemDefinition.EvaluateNormalPressureFunction(this->mCurrentTime);
1384 break;
1385 default:
1387 }
1388
1389 for (unsigned quad_index=0; quad_index<this->mpBoundaryQuadratureRule->GetNumQuadPoints(); quad_index++)
1390 {
1391 double wJ = jacobian_determinant * this->mpBoundaryQuadratureRule->GetWeight(quad_index);
1392
1393 // Get the quadrature point on this surface element (in canonical space) - so eg, for a 2D problem,
1394 // the quad point is in 1D space
1395 const ChastePoint<DIM-1>& quadrature_point = this->mpBoundaryQuadratureRule->rGetQuadPoint(quad_index);
1396 QuadraticBasisFunction<DIM-1>::ComputeBasisFunctions(quadrature_point, quad_phi_surf_element);
1397
1398 // We will need the xi coordinates of this quad point in the volume element. We could do this by figuring
1399 // out how the nodes of the surface element are ordered in the list of nodes in the volume element,
1400 // however it is less fiddly to compute directly. Firstly, compute the corresponding physical location
1401 // of the quad point, by interpolating
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++)
1404 {
1405 X += quad_phi_surf_element(node_index)*rBoundaryElement.GetNode(node_index)->rGetLocation();
1406 }
1407
1408
1409 // Now compute the xi coordinates of the quad point in the volume element
1410 c_vector<double,DIM+1> weight = p_containing_vol_element->CalculateInterpolationWeights(X);
1411 c_vector<double,DIM> xi;
1412 for (unsigned i=0; i<DIM; i++)
1413 {
1414 xi(i) = weight(i+1); // Note, in 2d say, weights = [1-xi(0)-xi(1), xi(0), xi(1)]
1415 }
1416
1417 // Check one of the weights was zero, as the quad point is on the boundary of the volume element
1418 if (DIM == 2)
1419 {
1420 assert( DIM!=2 || (fabs(weight(0))<1e-6) || (fabs(weight(1))<1e-6) || (fabs(weight(2))<1e-6) );
1421 }
1422 else
1423 {
1424 assert( DIM!=3 || (fabs(weight(0))<1e-6) || (fabs(weight(1))<1e-6) || (fabs(weight(2))<1e-6) || (fabs(weight(3))<1e-6) ); // LCOV_EXCL_LINE
1425 }
1426
1427 // Now we can compute the grad_phi and then interpolate F
1428 QuadraticBasisFunction<DIM>::ComputeTransformedBasisFunctionDerivatives(xi, inverse_jacobian_vol_element, grad_quad_phi_vol_element);
1429
1430 F = identity_matrix<double>(DIM,DIM);
1431 for (unsigned node_index=0; node_index<NUM_NODES_PER_ELEMENT; node_index++)
1432 {
1433 for (unsigned i=0; i<DIM; i++)
1434 {
1435 for (unsigned M=0; M<DIM; M++)
1436 {
1437 F(i,M) += grad_quad_phi_vol_element(M,node_index)*element_current_displacements(i,node_index);
1438 }
1439 }
1440 }
1441
1442 double detF = Determinant(F);
1443 invF = Inverse(F);
1444
1445 if (assembleResidual)
1446 {
1447 c_vector<double,DIM> traction = detF*normal_pressure*prod(trans(invF),normal);
1448
1449 // assemble
1450 for (unsigned index=0; index<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index++)
1451 {
1452 unsigned spatial_dim = index%DIM;
1453 unsigned node_index = (index-spatial_dim)/DIM;
1454
1455 assert(node_index < NUM_NODES_PER_BOUNDARY_ELEMENT);
1456
1457 rBelem(index) -= traction(spatial_dim)
1458 * quad_phi_surf_element(node_index)
1459 * wJ;
1460 }
1461 }
1462
1463 // Sometimes we don't include the analytic jacobian for this integral
1464 // in the jacobian that is assembled - the ShouldAssembleMatrixTermForPressureOnDeformedBc()
1465 // bit below - see the documentation for this method to see why.
1466 if (assembleJacobian && ShouldAssembleMatrixTermForPressureOnDeformedBc())
1467 {
1468 for (unsigned II=0; II<NUM_NODES_PER_BOUNDARY_ELEMENT; II++)
1469 {
1470 for (unsigned N=0; N<DIM; N++)
1471 {
1472 grad_quad_phi_surf_element(N,II) = grad_quad_phi_vol_element(N,surf_to_vol_map[II]);
1473 }
1474 }
1475
1477 for (unsigned N=0; N<DIM; N++)
1478 {
1479 for (unsigned e=0; e<DIM; e++)
1480 {
1481 for (unsigned M=0; M<DIM; M++)
1482 {
1483 for (unsigned d=0; d<DIM; d++)
1484 {
1485 tensor1(N,e,M,d) = invF(N,e)*invF(M,d) - invF(M,e)*invF(N,d);
1486 }
1487 }
1488 }
1489 }
1490
1491 // tensor2(II,e,M,d) = tensor1(N,e,M,d)*grad_quad_phi_surf_element(N,II)
1493 tensor2.template SetAsContractionOnFirstDimension<DIM>( trans(grad_quad_phi_surf_element), tensor1);
1494
1495 // tensor3 is really a third-order tensor
1496 // tensor3(II,e,0,d) = tensor2(II,e,M,d)*normal(M)
1498 tensor3.template SetAsContractionOnThirdDimension<DIM>( normal_as_mat, tensor2);
1499
1500 for (unsigned index1=0; index1<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index1++)
1501 {
1502 unsigned spatial_dim1 = index1%DIM;
1503 unsigned node_index1 = (index1-spatial_dim1)/DIM;
1504
1505 for (unsigned index2=0; index2<NUM_NODES_PER_BOUNDARY_ELEMENT*DIM; index2++)
1506 {
1507 unsigned spatial_dim2 = index2%DIM;
1508 unsigned node_index2 = (index2-spatial_dim2)/DIM;
1509
1510 rAelem(index1,index2) -= normal_pressure
1511 * detF
1512 * tensor3(node_index2,spatial_dim2,0,spatial_dim1)
1513 * quad_phi_surf_element(node_index1)
1514 * wJ;
1515 }
1516 }
1517 }
1518 }
1519}
1520
1521template<unsigned DIM>
1523{
1524 // Check the problem definition is set up correctly (and fully).
1525 mrProblemDefinition.Validate();
1526
1527 // If the problem includes specified pressures on deformed surfaces (as opposed
1528 // to specified tractions), the code needs to compute normals, and they need
1529 // to be consistently all facing outward (or all facing inward). Check the undeformed
1530 // mesh boundary elements has nodes that are ordered so that all normals are
1531 // outward-facing
1532 if (mrProblemDefinition.GetTractionBoundaryConditionType()==PRESSURE_ON_DEFORMED && mCheckedOutwardNormals==false)
1533 {
1534 this->mrQuadMesh.CheckOutwardNormals();
1535 mCheckedOutwardNormals = true;
1536 }
1537
1538 // Write the initial solution
1539 this->WriteCurrentSpatialSolution("initial", "nodes");
1540
1541 if (mUseSnesSolver)
1542 {
1543 SolveSnes();
1544 }
1545 else
1546 {
1547 SolveNonSnes(tol);
1548 }
1549
1550 // Remove pressure dummy values (P=0 at internal nodes, which should have been
1551 // been the result of the solve above), by linear interpolating from vertices of
1552 // edges to the internal node of the edge
1553 if (this->mCompressibilityType==INCOMPRESSIBLE)
1554 {
1555 this->RemovePressureDummyValuesThroughLinearInterpolation();
1556 }
1557
1558 // Write the final solution
1559 this->WriteCurrentSpatialSolution("solution", "nodes");
1560}
1561
1565template<unsigned DIM>
1567{
1568 // Four alternatives
1569 // (a) Petsc direct solve
1570 // Otherwise iterative solve with:
1571 // (b) Incompressible: GMRES with ILU preconditioner (or bjacobi=ILU on each process) [default]. Very poor on large problems.
1572 // (c) Incompressible: GMRES with AMG preconditioner. Uncomment #define MECH_USE_HYPRE above. Requires Petsc3 with HYPRE installed.
1573 // (d) Compressible: CG with ICC
1574
1575 PC pc;
1576 KSPGetPC(solver, &pc);
1577
1578 if (mPetscDirectSolve)
1579 {
1580 if (this->mCompressibilityType==COMPRESSIBLE)
1581 {
1582 KSPSetType(solver,KSPPREONLY);
1583
1584 }
1585 PCSetType(pc, PCLU);
1586
1587 // See #2057
1588 // PCFactorSetMatSolverPackage(pc,"mumps");
1589 }
1590 else
1591 {
1592 if (this->mCompressibilityType==COMPRESSIBLE)
1593 {
1594 KSPSetType(solver,KSPCG);
1596 {
1597 PCSetType(pc, PCICC);
1598 //Note that PetscOptionsSetValue is dangerous because we can't easily do
1599 //regression testing. If a name changes, then the behaviour of the code changes
1600 //because it won't recognise the old name. However, it won't fail to compile/run.
1601 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 1) //PETSc 3.1 or later
1602 PetscTools::SetOption("-pc_factor_shift_type", "positive_definite");
1603 #else
1604 PetscTools::SetOption("-pc_factor_shift_positive_definite", "");
1605 #endif
1606 }
1607 else
1608 {
1609 PCSetType(pc, PCBJACOBI);
1610 }
1611 }
1612 else
1613 {
1614 unsigned num_restarts = 100;
1615 KSPSetType(solver,KSPGMRES);
1616 KSPGMRESSetRestart(solver,num_restarts);
1617
1618 #ifndef MECH_USE_HYPRE
1619 PCSetType(pc, PCBJACOBI); // BJACOBI = ILU on each block (block = part of matrix on each process)
1620 #else
1622 // Speed up linear solve time massively for larger simulations (in fact GMRES may stagnate without
1623 // this for larger problems), by using a AMG preconditioner -- needs HYPRE installed
1625 PetscTools::SetOption("-pc_hypre_type", "boomeramg");
1626 // PetscTools::SetOption("-pc_hypre_boomeramg_max_iter", "1");
1627 // PetscTools::SetOption("-pc_hypre_boomeramg_strong_threshold", "0.0");
1628
1629 PCSetType(pc, PCHYPRE);
1630 #if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >=2) //PETSc 3.2 or later
1631 KSPSetPCSide(solver, PC_RIGHT);
1632 #else
1633 KSPSetPreconditionerSide(solver, PC_RIGHT);
1634 #endif
1635
1636 // other possible preconditioners..
1637 //PCBlockDiagonalMechanics* p_custom_pc = new PCBlockDiagonalMechanics(solver, this->mPreconditionMatrix, mBlock1Size, mBlock2Size);
1638 //PCLDUFactorisationMechanics* p_custom_pc = new PCLDUFactorisationMechanics(solver, this->mPreconditionMatrix, mBlock1Size, mBlock2Size);
1639 //remember to delete memory..
1640 //KSPSetPreconditionerSide(solver, PC_RIGHT);
1641 #endif
1642 }
1643 }
1644}
1645
1647// The code for the non-SNES solver - maybe remove all this
1648// as SNES solver appears better
1650
1651template<unsigned DIM>
1653{
1654 if (!allowException)
1655 {
1656 // Assemble the residual
1657 AssembleSystem(true, false);
1658 }
1659 else
1660 {
1661 try
1662 {
1663 // Try to assemble the residual using this current solution
1664 AssembleSystem(true, false);
1665 }
1666 catch(Exception&)
1667 {
1668 // If fail (because e.g. ODEs fail to solve, or strains are too large for material law), return infinity
1669 return DBL_MAX;
1670 }
1671 }
1672
1673 // Return the scaled norm of the residual
1674 return CalculateResidualNorm();
1675}
1676
1677template<unsigned DIM>
1679{
1680 double norm;
1681
1682 //\todo Change to NORM_1 and remove the division by mNumDofs...
1683 VecNorm(this->mResidualVector, NORM_2, &norm);
1684 return norm/this->mNumDofs;
1685}
1686
1687template<unsigned DIM>
1690 double a,
1691 std::vector<double>& rZ)
1692{
1693 assert(rX.size()==rY.GetSize());
1694 assert(rY.GetSize()==rZ.size());
1695 for (unsigned i=0; i<rX.size(); i++)
1696 {
1697 rZ[i] = rX[i] + a*rY[i];
1698 }
1699}
1700
1701template<unsigned DIM>
1703{
1704 if (this->mVerbose)
1705 {
1706 Timer::Reset();
1707 }
1708
1710 // Assemble Jacobian (and preconditioner)
1712 MechanicsEventHandler::BeginEvent(MechanicsEventHandler::ASSEMBLE);
1713 AssembleSystem(true, true);
1714 MechanicsEventHandler::EndEvent(MechanicsEventHandler::ASSEMBLE);
1715 if (this->mVerbose)
1716 {
1717 Timer::PrintAndReset("AssembleSystem");
1718 }
1719
1721 // Solve the linear system.
1723 MechanicsEventHandler::BeginEvent(MechanicsEventHandler::SOLVE);
1724
1725 Vec solution;
1726 VecDuplicate(this->mResidualVector,&solution);
1727
1728 KSP solver;
1729 KSPCreate(PETSC_COMM_WORLD,&solver);
1730
1731#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
1732 KSPSetOperators(solver, mrJacobianMatrix, this->mPreconditionMatrix);
1733#else
1734 KSPSetOperators(solver, mrJacobianMatrix, this->mPreconditionMatrix, DIFFERENT_NONZERO_PATTERN /*in precond between successive solves*/);
1735#endif
1736
1737 // Set the type of KSP solver (CG, GMRES etc) and preconditioner (ILU, HYPRE, etc)
1738 SetKspSolverAndPcType(solver);
1739
1740 //PetscTools::SetOption("-ksp_monitor","");
1741 //PetscTools::SetOption("-ksp_norm_type","natural");
1742
1743 KSPSetFromOptions(solver);
1744 KSPSetUp(solver);
1745
1746
1747 // Set the linear system absolute tolerance.
1748 // This is either the user provided value, or set to
1749 // max {1e-6 * initial_residual, 1e-12}
1750 if (mKspAbsoluteTol < 0)
1751 {
1752 Vec temp;
1753 VecDuplicate(this->mResidualVector, &temp);
1754 Vec temp2;
1755 VecDuplicate(this->mResidualVector, &temp2);
1756 Vec linsys_residual;
1757 VecDuplicate(this->mResidualVector, &linsys_residual);
1758
1759 KSPInitialResidual(solver, solution, temp, temp2, linsys_residual, this->mLinearSystemRhsVector);
1760 double initial_resid_norm;
1761 VecNorm(linsys_residual, NORM_2, &initial_resid_norm);
1762
1763 PetscTools::Destroy(temp);
1764 PetscTools::Destroy(temp2);
1765 PetscTools::Destroy(linsys_residual);
1766
1767 double ksp_rel_tol = 1e-6;
1768 double absolute_tol = ksp_rel_tol * initial_resid_norm;
1769 if (absolute_tol < 1e-12)
1770 {
1771 absolute_tol = 1e-12;
1772 }
1773 KSPSetTolerances(solver, 1e-16, absolute_tol, PETSC_DEFAULT, 1000 /* max iters */); // Note: some machines - max iters seems to be 1000 whatever we give here
1774 }
1775 else
1776 {
1777 KSPSetTolerances(solver, 1e-16, mKspAbsoluteTol, PETSC_DEFAULT, 1000 /* max iters */); // Note: some machines - max iters seems to be 1000 whatever we give here
1778 }
1779
1780 if (this->mVerbose)
1781 {
1782 Timer::PrintAndReset("KSP Setup");
1783 }
1784
1785 KSPSolve(solver,this->mLinearSystemRhsVector,solution);
1786
1787// ///// For printing matrix when debugging
1788// OutputFileHandler handler("TEMP",false);
1789// std::stringstream ss;
1790// static unsigned counter = 0;
1791// ss << "all_" << counter++ << ".txt";
1792// out_stream p_file = handler.OpenOutputFile(ss.str());
1793// *p_file << std::setprecision(10);
1794// for (unsigned i=0; i<this->mNumDofs; i++)
1795// {
1796// for (unsigned j=0; j<this->mNumDofs; j++)
1797// {
1798// *p_file << PetscMatTools::GetElement(mrJacobianMatrix, i, j) << " ";
1799// }
1800// *p_file << PetscVecTools::GetElement(this->mLinearSystemRhsVector, i) << " ";
1801// *p_file << PetscVecTools::GetElement(solution, i) << "\n";
1802// }
1803// p_file->close();
1804
1805
1807 // Error checking for linear solve
1809
1810 // warn if ksp reports failure
1811 KSPConvergedReason reason;
1812 KSPGetConvergedReason(solver,&reason);
1813 if (reason == KSP_DIVERGED_ITS || reason == KSP_DIVERGED_BREAKDOWN)
1814 {
1815 // DIVERGED_ITS or DIVERGED_BREAKDOWN just means it didn't converge in the given maximum number of iterations,
1816 // or similar, which is potentially not a problem, as the nonlinear solver may (and often will) still converge.
1817 // Just warn once.
1818 // (Very difficult to cover in normal tests, requires relative and absolute ksp tols to be very small, there
1819 // is no interface for setting both of these. Could be covered by setting up a problem the solver
1820 // finds difficult to solve, but this would be overkill.)
1821 // LCOV_EXCL_START
1822 WARN_ONCE_ONLY("Linear solve (within a mechanics solve) didn't converge, but this may not stop nonlinear solve converging")
1823 // LCOV_EXCL_STOP
1824 }
1825 else
1826 {
1827 // Throw an exception if the solver failed for any reason other than DIVERGED_ITS.
1828 // This is not covered as would be difficult to cover - requires a bad matrix to
1829 // assembled, for example.
1830 // LCOV_EXCL_START
1831 KSPEXCEPT(reason);
1832 // LCOV_EXCL_STOP
1833 }
1834
1835 // quit if no ksp iterations were done
1836 int num_iters;
1837 KSPGetIterationNumber(solver, &num_iters);
1838 if (num_iters==0)
1839 {
1840 PetscTools::Destroy(solution);
1841 KSPDestroy(PETSC_DESTROY_PARAM(solver));
1842 EXCEPTION("KSP Absolute tolerance was too high, linear system wasn't solved - there will be no decrease in Newton residual. Decrease KspAbsoluteTolerance");
1843 }
1844
1845
1846 if (this->mVerbose)
1847 {
1848 Timer::PrintAndReset("KSP Solve");
1849 std::cout << "[" << PetscTools::GetMyRank() << "]: Num iterations = " << num_iters << "\n" << std::flush;
1850 }
1851
1852 MechanicsEventHandler::EndEvent(MechanicsEventHandler::SOLVE);
1853
1855 // Update the solution
1856 // Newton method: sol = sol - update, where update=Jac^{-1}*residual
1857 // Newton with damping: sol = sol - s*update, where s is chosen
1858 // such that |residual(sol)| is minimised. Damping is important to
1859 // avoid initial divergence.
1860 //
1861 // Normally, finding the best s from say 0.05,0.1,0.2,..,1.0 is cheap,
1862 // but this is not the case in cardiac electromechanics calculations.
1863 // Therefore, we initially check s=1 (expected to be the best most of the
1864 // time, then s=0.9. If the norm of the residual increases, we assume
1865 // s=1 is the best. Otherwise, check s=0.8 to see if s=0.9 is a local min.
1867 MechanicsEventHandler::BeginEvent(MechanicsEventHandler::UPDATE);
1868 double new_norm_resid = UpdateSolutionUsingLineSearch(solution);
1869 MechanicsEventHandler::EndEvent(MechanicsEventHandler::UPDATE);
1870
1871 PetscTools::Destroy(solution);
1872 KSPDestroy(PETSC_DESTROY_PARAM(solver));
1873
1874 return new_norm_resid;
1875}
1876
1877template<unsigned DIM>
1879{
1880 if (this->mVerbose)
1881 {
1882 std::cout << "\tTesting s = " << s << ", |f| = " << residNorm << "\n" << std::flush;
1883 }
1884}
1885
1886template<unsigned DIM>
1888{
1889 double initial_norm_resid = CalculateResidualNorm();
1890 if (this->mVerbose)
1891 {
1892 std::cout << "\tInitial |f| [corresponding to s=0] is " << initial_norm_resid << "\n" << std::flush;
1893 }
1894
1895 ReplicatableVector update(solution);
1896
1897 std::vector<double> old_solution = this->mCurrentSolution;
1898
1899 std::vector<double> damping_values; // = {1.0, 0.9, .., 0.2, 0.1, 0.05} ie size 11
1900 for (unsigned i=10; i>=1; i--)
1901 {
1902 damping_values.push_back((double)i/10.0);
1903 }
1904 damping_values.push_back(0.05);
1905 assert(damping_values.size()==11);
1906
1908 // let mCurrentSolution = old_solution - damping_val[0]*update; and compute residual
1909 unsigned index = 0;
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);
1913
1915 // let mCurrentSolution = old_solution - damping_val[1]*update; and compute residual
1916 index = 1;
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);
1920
1921 index = 2;
1922 // While f(s_next) < f(s_current), [f = residnorm], keep trying new damping values,
1923 // ie exit thus loop when next norm of the residual first increases
1924 while ( (next_resid_norm==DBL_MAX) // the residual is returned as infinity if the deformation is so large to cause exceptions in the material law/EM contraction model
1925 || ( (next_resid_norm < current_resid_norm) && (index<damping_values.size()) ) )
1926 {
1927 current_resid_norm = next_resid_norm;
1928
1929 // let mCurrentSolution = old_solution - damping_val*update; and compute residual
1930 VectorSum(old_solution, update, -damping_values[index], this->mCurrentSolution);
1931 next_resid_norm = ComputeResidualAndGetNorm(true);
1932 PrintLineSearchResult(damping_values[index], next_resid_norm);
1933
1934 index++;
1935 }
1936
1937 unsigned best_index;
1938
1939 if (index==damping_values.size() && (next_resid_norm < current_resid_norm))
1940 {
1941 // Difficult to come up with large forces/tractions such that it had to
1942 // test right down to s=0.05, but overall doesn't fail.
1943 // The possible damping values have been manually temporarily altered to
1944 // get this code to be called, it appears to work correctly. Even if it
1945 // didn't tests wouldn't fail, they would just be v. slightly less efficient.
1946 // LCOV_EXCL_START
1947 // if we exited because we got to the end of the possible damping values, the
1948 // best one was the last one (excl the final index++ at the end)
1949 current_resid_norm = next_resid_norm;
1950 best_index = index-1;
1951 // LCOV_EXCL_STOP
1952 }
1953 else
1954 {
1955 // else the best one must have been the second last one (excl the final index++ at the end)
1956 // (as we would have exited when the resid norm first increased)
1957 best_index = index-2;
1958 }
1959
1960 // See documentation for SetTakeFullFirstNewtonStep()
1961 bool full_first_step = mTakeFullFirstNewtonStep && mFirstStep;
1962
1963
1964 // Check out best was better than the original residual-norm
1965 if (initial_norm_resid < current_resid_norm && !full_first_step)
1966 {
1967 // LCOV_EXCL_START
1968 EXCEPTION("Residual does not appear to decrease in newton direction, quitting");
1969 // LCOV_EXCL_STOP
1970 }
1971
1972 // See documentation for SetTakeFullFirstNewtonStep()
1973 if (full_first_step)
1974 {
1975 if (this->mVerbose)
1976 {
1977 std::cout << "\tTaking full first Newton step...\n";
1978 }
1979 best_index = 0;
1980 }
1981
1982 if (this->mVerbose)
1983 {
1984 std::cout << "\tChoosing s = " << damping_values[best_index] << "\n" << std::flush;
1985 }
1986
1987
1989//
1990// double l_inf_disp = 0.0;
1991// double l_inf_pressure = 0.0;
1992//
1993// if (this->mCompressibilityType==INCOMPRESSIBLE)
1994// {
1995// for (unsigned i=0; i<this->mrQuadMesh.GetNumNodes(); i++)
1996// {
1997// for (unsigned j=0; j<DIM; j++)
1998// {
1999// double value = update[(DIM+1)*i + j]*damping_values[best_index];
2000// l_inf_disp = std::max(l_inf_disp, fabs(value));
2001// }
2002// l_inf_pressure = std::max(l_inf_pressure, fabs(update[(DIM+1)*i + DIM]*damping_values[best_index]));
2003// }
2004// std::cout << "l_inf_disp, l_inf_pressure = " << l_inf_disp << " " << l_inf_pressure << "\n";
2005// }
2006// else
2007// {
2008// for (unsigned i=0; i<this->mrQuadMesh.GetNumNodes(); i++)
2009// {
2010// for (unsigned j=0; j<DIM; j++)
2011// {
2012// double value = update[DIM*i + j]*damping_values[best_index];
2013// l_inf_disp = std::max(l_inf_disp, fabs(value));
2014// }
2015// }
2016// std::cout << "l_inf_disp = " << l_inf_disp << "\n";
2017// }
2018
2019 VectorSum(old_solution, update, -damping_values[best_index], this->mCurrentSolution);
2020
2021 mLastDampingValue = damping_values[best_index];
2022 return current_resid_norm;
2023}
2024
2025template<unsigned DIM>
2026void AbstractNonlinearElasticitySolver<DIM>::PostNewtonStep(unsigned counter, double normResidual)
2027{
2028}
2029
2030template<unsigned DIM>
2032{
2033 mLastDampingValue = 0;
2034
2035 if (mWriteOutputEachNewtonIteration)
2036 {
2037 this->WriteCurrentSpatialSolution("newton_iteration", "nodes", 0);
2038 }
2039
2040 // Compute residual
2041 double norm_resid = ComputeResidualAndGetNorm(false);
2042 if (this->mVerbose)
2043 {
2044 std::cout << "\nNorm of residual is " << norm_resid << "\n";
2045 }
2046
2047 mNumNewtonIterations = 0;
2048 unsigned iteration_number = 1;
2049
2050 if (tol < 0) // i.e. if wasn't passed in as a parameter
2051 {
2052 tol = NEWTON_REL_TOL*norm_resid;
2053
2054 // LCOV_EXCL_START // not going to have tests in cts for everything
2055 if (tol > MAX_NEWTON_ABS_TOL)
2056 {
2057 tol = MAX_NEWTON_ABS_TOL;
2058 }
2059 if (tol < MIN_NEWTON_ABS_TOL)
2060 {
2061 tol = MIN_NEWTON_ABS_TOL;
2062 }
2063 // LCOV_EXCL_STOP
2064 }
2065
2066 if (this->mVerbose)
2067 {
2068 std::cout << "Solving with tolerance " << tol << "\n";
2069 }
2070
2071 while (norm_resid > tol)
2072 {
2073 if (this->mVerbose)
2074 {
2075 std::cout << "\n-------------------\n"
2076 << "Newton iteration " << iteration_number
2077 << ":\n-------------------\n";
2078 }
2079
2080 // take newton step (and get returned residual)
2081 mFirstStep = (iteration_number==1);
2082 norm_resid = TakeNewtonStep();
2083
2084 if (this->mVerbose)
2085 {
2086 std::cout << "Norm of residual is " << norm_resid << "\n";
2087 }
2088
2089 if (mWriteOutputEachNewtonIteration)
2090 {
2091 this->WriteCurrentSpatialSolution("newton_iteration", "nodes", iteration_number);
2092 }
2093
2094 mNumNewtonIterations = iteration_number;
2095
2096 PostNewtonStep(iteration_number,norm_resid);
2097
2098 iteration_number++;
2099 if (iteration_number==20)
2100 {
2101 // LCOV_EXCL_START
2102 EXCEPTION("Not converged after 20 newton iterations, quitting");
2103 // LCOV_EXCL_STOP
2104 }
2105 }
2106
2107 if (norm_resid > tol)
2108 {
2109 // LCOV_EXCL_START
2110 EXCEPTION("Failed to converge");
2111 // LCOV_EXCL_STOP
2112 }
2113}
2114
2115template<unsigned DIM>
2117{
2118 return mNumNewtonIterations;
2119}
2120
2122// SNES version of the nonlinear solver
2124
2125template<unsigned DIM>
2127{
2128 // Set up solution guess for residuals
2129 Vec initial_guess;
2130 VecDuplicate(this->mResidualVector, &initial_guess);
2131 double* p_initial_guess;
2132 VecGetArray(initial_guess, &p_initial_guess);
2133 int lo, hi;
2134 VecGetOwnershipRange(initial_guess, &lo, &hi);
2135 for (int global_index=lo; global_index<hi; global_index++)
2136 {
2137 int local_index = global_index - lo;
2138 p_initial_guess[local_index] = this->mCurrentSolution[global_index];
2139 }
2140 VecRestoreArray(initial_guess, &p_initial_guess);
2141 PetscVecTools::Finalise(initial_guess);
2142
2143 Vec snes_residual_vec;
2144 VecDuplicate(this->mResidualVector, &snes_residual_vec);
2145
2146 SNES snes;
2147
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) //PETSc 3.4 or later
2152 SNESSetType(snes, SNESNEWTONLS);
2153#else
2154 SNESSetType(snes, SNESLS);
2155#endif
2156 SNESSetTolerances(snes, 1e-5, 1e-5, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT);
2157
2158#if (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR == 3) //PETSc 3.3
2159 SNESLineSearch linesearch;
2160 SNESGetSNESLineSearch(snes, &linesearch);
2161 // Use 'critical point' line search algorithm. This was changed from 'backtracking'; see #2916
2162 SNESLineSearchSetType(linesearch, "cp");
2163#elif (PETSC_VERSION_MAJOR == 3 && PETSC_VERSION_MINOR >= 4) //PETSc 3.4 or later
2164 SNESLineSearch linesearch;
2165 SNESGetLineSearch(snes, &linesearch);
2166 // Use 'critical point' line search algorithm. This was changed from 'backtracking'; see #2916
2167 SNESLineSearchSetType(linesearch, "cp");
2168#endif
2169
2170 SNESSetMaxLinearSolveFailures(snes,100);
2171
2172 KSP ksp;
2173 SNESGetKSP(snes, &ksp);
2174
2175 KSPSetTolerances(ksp, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1000 /* max iters */); // Note: some machines - max iters seems to be 1000 whatever we give here
2176
2177 // Set the type of KSP solver (CG, GMRES etc) and preconditioner (ILU, HYPRE, etc)
2178 SetKspSolverAndPcType(ksp);
2179
2180 if (this->mVerbose)
2181 {
2182 PetscTools::SetOption("-snes_monitor","");
2183 }
2184 SNESSetFromOptions(snes);
2185
2186#if (PETSC_VERSION_MAJOR == 2 && PETSC_VERSION_MINOR == 2) //PETSc 2.2
2187 PetscErrorCode err = SNESSolve(snes, initial_guess);
2188#else
2189 PetscErrorCode err = SNESSolve(snes, CHASTE_PETSC_NULLPTR, initial_guess);
2190#endif
2191
2192 SNESConvergedReason reason;
2193 SNESGetConvergedReason(snes,&reason);
2194
2195// LCOV_EXCL_START
2196 if (err != 0)
2197 {
2198 std::stringstream err_stream;
2199 err_stream << err;
2200 PetscTools::Destroy(initial_guess);
2201 PetscTools::Destroy(snes_residual_vec);
2202 SNESDestroy(PETSC_DESTROY_PARAM(snes));
2203 EXCEPTION("Nonlinear Solver failed. PETSc error code: "+err_stream.str()+" .");
2204 }
2205
2206 if (reason < 0)
2207 {
2208 std::stringstream reason_stream;
2209 reason_stream << reason;
2210 PetscTools::Destroy(initial_guess);
2211 PetscTools::Destroy(snes_residual_vec);
2212 SNESDestroy(PETSC_DESTROY_PARAM(snes));
2213 EXCEPTION("Nonlinear Solver did not converge. PETSc reason code: "+reason_stream.str()+" .");
2214 }
2215// LCOV_EXCL_STOP
2216
2217 PetscInt num_iters;
2218 SNESGetIterationNumber(snes,&num_iters);
2219 mNumNewtonIterations = num_iters;
2220
2221 PetscTools::Destroy(initial_guess);
2222 PetscTools::Destroy(snes_residual_vec);
2223 SNESDestroy(PETSC_DESTROY_PARAM(snes));
2224}
2225
2226template<unsigned DIM>
2228{
2229 // Note: AssembleSystem() assumes the current solution is in this->mCurrentSolution and assembles
2230 // this->mResiduaVector and/or this->mrJacobianMatrix. Since PETSc wants us to use the input
2231 // currentGuess, and write the output to residualVector, we have to copy do some copies below.
2232
2233 ReplicatableVector guess_repl(currentGuess);
2234 for (unsigned i=0; i<guess_repl.GetSize(); i++)
2235 {
2236 this->mCurrentSolution[i] = guess_repl[i];
2237 }
2238 AssembleSystem(true,false);
2239 VecCopy(this->mResidualVector, residualVector);
2240}
2241
2242template<unsigned DIM>
2243void AbstractNonlinearElasticitySolver<DIM>::ComputeJacobian(Vec currentGuess, Mat* pJacobian, Mat* pPreconditioner)
2244{
2245 // Note: AssembleSystem() assumes the current solution is in this->mCurrentSolution and assembles
2246 // this->mResiduaVector and/or this->mrJacobianMatrix.
2247 // We need to copy the input currentGuess into the local mCurrentGuess.
2248 // We don't have to copy mrJacobianMatrix to pJacobian, which would be expensive, as they will
2249 // point to the same memory.
2250
2251 // check Petsc data corresponds to internal Mats
2252 assert(mrJacobianMatrix==*pJacobian);
2253 assert(this->mPreconditionMatrix==*pPreconditioner);
2254
2255 MechanicsEventHandler::BeginEvent(MechanicsEventHandler::ASSEMBLE);
2256 ReplicatableVector guess_repl(currentGuess);
2257 for (unsigned i=0; i<guess_repl.GetSize(); i++)
2258 {
2259 this->mCurrentSolution[i] = guess_repl[i];
2260 }
2261
2262 AssembleSystem(false,true);
2263 MechanicsEventHandler::EndEvent(MechanicsEventHandler::ASSEMBLE);
2264}
2265
2266template<unsigned DIM>
2267PetscErrorCode AbstractNonlinearElasticitySolver_ComputeResidual(SNES snes,
2268 Vec currentGuess,
2269 Vec residualVector,
2270 void* pContext)
2271{
2272 // Extract the solver from the void*
2274 p_solver->ComputeResidual(currentGuess, residualVector);
2275 return 0;
2276}
2277
2278template<unsigned DIM>
2279#if ((PETSC_VERSION_MAJOR==3) && (PETSC_VERSION_MINOR>=5))
2280 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
2281 Vec currentGuess,
2282 Mat globalJacobian,
2283 Mat preconditioner,
2284 void* pContext)
2285 {
2286 // Extract the solver from the void*
2288 p_solver->ComputeJacobian(currentGuess, &globalJacobian, &preconditioner);
2289 return 0;
2290 }
2291#else
2292 PetscErrorCode AbstractNonlinearElasticitySolver_ComputeJacobian(SNES snes,
2293 Vec currentGuess,
2294 Mat* pGlobalJacobian,
2295 Mat* pPreconditioner,
2296 MatStructure* pMatStructure,
2297 void* pContext)
2298 {
2299 // Extract the solver from the void*
2301 p_solver->ComputeJacobian(currentGuess, pGlobalJacobian, pPreconditioner);
2302 return 0;
2303 }
2304#endif
2305
2306
2307// Constant setting definitions
2308
2309template<unsigned DIM>
2311
2312template<unsigned DIM>
2314
2315template<unsigned DIM>
2317
2318#endif /*ABSTRACTNONLINEARELASTICITYSOLVER_HPP_*/
#define EXCEPTION(message)
#define NEVER_REACHED
#define CHASTE_PETSC_NULLPTR
A macro to define PETSc null pointer based on the PETSc version.
#define PETSC_DESTROY_PARAM(x)
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
std::vector< c_vector< double, DIM > > & rGetSpatialSolution()
void SetIncludeActiveTension(bool includeActiveTension=true)
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)
std::vector< c_vector< double, DIM *(DIM+1)/2 > > mAverageStressesPerElement
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)
virtual void FinishAssembleSystem(bool assembleResidual, bool assembleLinearSystem)
void SetComputeAverageStressPerElementDuringSolve(bool setComputeAverageStressPerElement=true)
void ComputeResidual(Vec currentGuess, Vec residualVector)
FourthOrderTensor< DIM, DIM, DIM, DIM > dTdE
std::vector< c_vector< double, DIM > > & rGetDeformedPosition()
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)
void SetKspAbsoluteTolerance(double kspAbsoluteTolerance)
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)
SolidMechanicsProblemDefinition< DIM > & mrProblemDefinition
void WriteData(OutputFileHandler &rHandler, const std::string &rFileName)
c_vector< double, SPACE_DIM > CalculateNormal()
c_vector< double, DIM > & rGetLocation()
void WriteCmguiScript(std::string fieldBaseName="", std::string undeformedBaseName="")
void WriteDeformationPositions(std::vector< c_vector< double, DIM > > &rDeformedPositions, unsigned counter)
void WriteInitialMesh(std::string fileName="")
bool OptionExists(const std::string &rOption)
static CommandLineArguments * Instance()
c_vector< double, SPACE_DIM+1 > CalculateInterpolationWeights(const ChastePoint< SPACE_DIM > &rTestPoint)
Definition Element.cpp:224
static void SwitchWriteMode(Mat matrix)
static void Finalise(Mat matrix)
static void Destroy(Vec &rVec)
static bool IsSequential()
static void SetOption(const char *pOptionName, const char *pOptionValue)
static unsigned GetMyRank()
static void Finalise(Vec vector)
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)
Definition Timer.cpp:70
static void Reset()
Definition Timer.cpp:44