AbstractIsotropicIncompressibleMaterialLaw.cpp
00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029 #include "AbstractIsotropicIncompressibleMaterialLaw.hpp"
00030
00031 template<unsigned DIM>
00032 AbstractIsotropicIncompressibleMaterialLaw<DIM>::~AbstractIsotropicIncompressibleMaterialLaw()
00033 {
00034 }
00035
00036 template<unsigned DIM>
00037 void AbstractIsotropicIncompressibleMaterialLaw<DIM>::ComputeStressAndStressDerivative(
00038 c_matrix<double,DIM,DIM>& rC,
00039 c_matrix<double,DIM,DIM>& rInvC,
00040 double pressure,
00041 c_matrix<double,DIM,DIM>& rT,
00042 FourthOrderTensor<DIM,DIM,DIM,DIM>& rDTdE,
00043 bool computeDTdE)
00044 {
00045
00046
00047
00048
00049 #define COVERAGE_IGNORE
00050 assert((DIM==2) || (DIM==3));
00051 #undef COVERAGE_IGNORE
00052
00053 static c_matrix<double,DIM,DIM> identity = identity_matrix<double>(DIM);
00054
00055 double I1 = Trace(rC);
00056 double I2 = SecondInvariant(rC);
00057
00058 double w1 = Get_dW_dI1(I1, I2);
00059 double w2;
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 rT = 2*w1*identity - pressure*rInvC;
00070 if (DIM==3)
00071 {
00072 w2 = Get_dW_dI2(I1, I2);
00073 rT += 2*w2*(I1*identity - rC);
00074 }
00075
00076
00077
00078
00079
00080
00081
00082
00083
00084
00085
00086
00087
00088
00089
00090
00091
00092
00093
00094
00095
00096
00097
00098 if (computeDTdE)
00099 {
00100 double w11 = Get_d2W_dI1(I1,I2);
00101
00102 double w12;
00103 double w22;
00104
00105 if (DIM==3)
00106 {
00107 w22 = Get_d2W_dI2(I1, I2);
00108 w12 = Get_d2W_dI1I2(I1, I2);
00109 }
00110
00111 for (unsigned M=0; M<DIM; M++)
00112 {
00113 for (unsigned N=0; N<DIM; N++)
00114 {
00115 for (unsigned P=0; P<DIM; P++)
00116 {
00117 for (unsigned Q=0; Q<DIM; Q++)
00118 {
00119 rDTdE(M,N,P,Q) = 4 * w11 * (M==N) * (P==Q)
00120 + 2 * pressure * rInvC(M,P) * rInvC(Q,N);
00121
00122 if (DIM==3)
00123 {
00124 rDTdE(M,N,P,Q) += 4 * w22 * (I1*(M==N) - rC(M,N)) * (I1*(P==Q) - rC(P,Q))
00125 + 4 * w2 * ((M==N)*(P==Q) - (M==P)*(N==Q))
00126 + 4 * w12 * ((M==N)*(I1*(P==Q) - rC(P,Q)) + (P==Q)*(I1*(M==N) - rC(M,N)));
00127 }
00128 }
00129 }
00130 }
00131 }
00132 }
00133 }
00134
00135 template<>
00136 double AbstractIsotropicIncompressibleMaterialLaw<2>::GetZeroStrainPressure()
00137 {
00138 return 2*Get_dW_dI1(2,0);
00139 }
00140
00141 template<>
00142 double AbstractIsotropicIncompressibleMaterialLaw<3>::GetZeroStrainPressure()
00143 {
00144 return 2*Get_dW_dI1(3,3) + 4*Get_dW_dI2(3,3);
00145 }
00146
00148
00150
00151 template class AbstractIsotropicIncompressibleMaterialLaw<2>;
00152 template class AbstractIsotropicIncompressibleMaterialLaw<3>;