45 std::vector<std::vector<double> > a,
46 std::vector<std::vector<double> > b)
48 if (DIM!=2 && DIM !=3)
50 EXCEPTION(
"Can only have a 2 or 3d incompressible pole-zero law");
53 assert(k.size() == DIM);
54 assert(a.size() == DIM);
55 assert(b.size() == DIM);
57 for (
unsigned i = 0; i < DIM; i++)
59 assert(k[i].size() == DIM);
60 assert(a[i].size() == DIM);
61 assert(b[i].size() == DIM);
63 for (
unsigned j = 0; j < DIM; j++)
65 assert(k[i][j] == k[j][i]);
66 assert(a[i][j] == a[j][i]);
67 assert(b[i][j] == b[j][i]);
75 for (
unsigned M=0; M<DIM; M++)
77 for (
unsigned N=0; N<DIM; N++)
79 mIdentity(M,N) = M==N ? 1.0 : 0.0;
94 c_matrix<double,DIM,DIM>& rInvC,
96 c_matrix<double,DIM,DIM>& rT,
100 static c_matrix<double,DIM,DIM> C_transformed;
101 static c_matrix<double,DIM,DIM> invC_transformed;
109 this->ComputeTransformedDeformationTensor(rC, rInvC, C_transformed, invC_transformed);
113 c_matrix<double,DIM,DIM> E = 0.5*(C_transformed - mIdentity);
115 for (
unsigned M=0; M<DIM; M++)
117 for (
unsigned N=0; N<DIM; N++)
128 EXCEPTION(
"E_{MN} >= a_{MN} - strain unacceptably large for model");
135 - pressure*invC_transformed(M,N);
142 for (
unsigned M=0; M<DIM; M++)
144 for (
unsigned N=0; N<DIM; N++)
146 for (
unsigned P=0; P<DIM; P++)
148 for (
unsigned Q=0; Q<DIM; Q++)
150 rDTdE(M,N,P,Q) = 2 * pressure * invC_transformed(M,P) * invC_transformed(Q,N);
172 this->TransformStressAndStressDerivative(rT, rDTdE, computeDTdE);
void ComputeStressAndStressDerivative(c_matrix< double, DIM, DIM > &rC, c_matrix< double, DIM, DIM > &rInvC, double pressure, c_matrix< double, DIM, DIM > &rT, FourthOrderTensor< DIM, DIM, DIM, DIM > &rDTdE, bool computeDTdE)