OOFEM  2.4
OOFEM.org - Object Oriented Finite Element Solver
mooneyrivlin.C
Go to the documentation of this file.
1 /*
2  *
3  * ##### ##### ###### ###### ### ###
4  * ## ## ## ## ## ## ## ### ##
5  * ## ## ## ## #### #### ## # ##
6  * ## ## ## ## ## ## ## ##
7  * ## ## ## ## ## ## ## ##
8  * ##### ##### ## ###### ## ##
9  *
10  *
11  * OOFEM : Object Oriented Finite Element Code
12  *
13  * Copyright (C) 1993 - 2014 Borek Patzak
14  *
15  *
16  *
17  * Czech Technical University, Faculty of Civil Engineering,
18  * Department of Structural Mechanics, 166 29 Prague, Czech Republic
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * This program is distributed in the hope that it will be useful,
26  * but WITHOUT ANY WARRANTY; without even the implied warranty of
27  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
28  * Lesser General Public License for more details.
29  *
30  * You should have received a copy of the GNU Lesser General Public
31  * License along with this library; if not, write to the Free Software
32  * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
33  */
34 
35 #include "mooneyrivlin.h"
36 #include "floatmatrix.h"
37 #include "floatarray.h"
38 #include "classfactory.h"
39 #include "mathfem.h"
40 
41 namespace oofem {
42 REGISTER_Material(MooneyRivlinMaterial);
43 
44 MooneyRivlinMaterial :: MooneyRivlinMaterial(int n, Domain *d) : StructuralMaterial(n, d)
45 { }
46 
47 
48 void
49 MooneyRivlinMaterial :: giveFirstPKStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &vF, TimeStep *tStep)
50 // returns 9 components of the first piola kirchhoff stress corresponding to the given deformation gradinet
51 
52 {
53  double J, I1, I2, barI1, barI2;
54  FloatMatrix F, C, Cpow2, invFt, FC, invF;
55 
56  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
57 
58  //store deformation gradient into matrix
59  F.beMatrixForm(vF);
60  invF.beInverseOf(F);
61  invFt.beTranspositionOf(invF);
62  // compute jacobian
63  J = F.giveDeterminant();
64  // compute right Cauchy-Green tensor
65  C.beTProductOf(F, F);
66  // compute C*C
67  Cpow2.beProductOf(C, C);
68  // compute F*C
69  FC.beProductOf(F, C);
70  //compute first invariant of deviatoric part of C;
71  I1 = ( C.at(1, 1) + C.at(2, 2) + C.at(3, 3) );
72  I2 = 1. / 2. * ( I1 * I1 - Cpow2.at(1, 1) - Cpow2.at(2, 2) - Cpow2.at(3, 3) );
73  barI1 = I1 * pow(J, -2. / 3.);
74  barI2 = I2 * pow(J, -4. / 3.);
75 
76 
77  FloatMatrix P;
78  //first part of stress tensor : C1 * \frac{\partial \bar{I}_1}{\partial F_ij}
79  P.add(2 * C1 / pow(J, 2 / 3.), F);
80  P.add(-2. / 3. * C1 * barI1, invFt);
81  // second part of stress tensor : C2 * \frac{\partial \bar{I}_2}{\partial F_ij}
82  P.add( 2. * C2 * barI1 / pow(J, 2. / 3.), F );
83  P.add(-4. / 3. * C2 * barI2, invFt);
84  P.add(-2. * C2 / pow(J, 4. / 3.), FC );
85  // third part of stress tensor : K * \frac{\partial ln J }{F_ij}
86  P.add(K, invFt);
87 
88  answer.beVectorForm(P);
89 
90  // update gp
91  status->letTempFVectorBe(vF);
92  status->letTempPVectorBe(answer);
93 }
94 
95 void
96 MooneyRivlinMaterial :: giveFirstPKStressVector_PlaneStrain(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedvF, TimeStep *tStep)
97 // returns 4 components of the First PK Stress corresponding to the given deformation gradient
98 
99 {
100  double J, I1, I2, barI1, barI2;
101  FloatArray vF, vP;
102  FloatMatrix F, C, CC, invFt, FC, invF;
103 
104  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
105  vF = {
106  reducedvF.at(1), reducedvF.at(2), reducedvF.at(3), 0, 0, reducedvF.at(4), 0, 0, reducedvF.at(5)
107  };
108  //store deformation gradient into matrix
109  F.beMatrixForm(vF);
110  invF.beInverseOf(F);
111  invFt.beTranspositionOf(invF);
112  // compute jacobian and its log
113  J = F.giveDeterminant();
114  // compute right Cauchy-Green tensor
115  C.beTProductOf(F, F);
116  // compute C*C
117  CC.beProductOf(C, C);
118  // compute F*C
119  FC.beProductOf(F, C);
120  //compute first invariant of deviatoric part of C;
121  I1 = ( C.at(1, 1) + C.at(2, 2) + C.at(3, 3) );
122  I2 = 1. / 2. * ( I1 * I1 - CC.at(1, 1) - CC.at(2, 2) - CC.at(3, 3) );
123  barI1 = I1 * pow(J, -2. / 3.);
124  barI2 = I2 * pow(J, -4. / 3.);
125 
126 
127  FloatMatrix P;
128  //first part of stress tensor : C1 * \frac{\partial \bar{I}_1}{\partial F_ij}
129  P.add( 2 * C1 / pow(J, 2 / 3.), F );
130  P.add(-2. / 3. * C1 * barI1, invFt);
131  // second part of stress tensor : C2 * \frac{\partial \bar{I}_2}{\partial F_ij}
132  P.add( 2. * C2 * barI1 / pow(J, 2. / 3.), F );
133  P.add(-4. / 3. * C2 * barI2, invFt);
134  P.add(-2. * C2 / pow(J, 4. / 3.), FC);
135  // third part of stress tensor : K * \frac{\partial ln J }{F_ij}
136  P.add(K, invFt);
137  vP.beVectorForm(P);
138  StructuralMaterial :: giveReducedVectorForm(answer, vP, _PlaneStrain);
139  answer.beVectorForm(P);
140 
141  // update gp
142  status->letTempFVectorBe(reducedvF);
143  status->letTempPVectorBe(answer);
144 }
145 
146 
147 void
148 MooneyRivlinMaterial :: give3dMaterialStiffnessMatrix_dPdF(FloatMatrix &answer,
149  MatResponseMode mode,
150  GaussPoint *gp, TimeStep *tStep)
151 // returns the 9x9 tangent stiffness matrix - dP/dF
152 {
153  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
154 
155  double I1, I2, J, lnJ;
156  FloatMatrix C, CC, invF, FC, B;
157 
158 
159  FloatArray vF = status->giveTempFVector();
160  FloatMatrix F = { { vF.at(1), vF.at(6), vF.at(5) }, { vF.at(9), vF.at(2), vF.at(4) }, { vF.at(8), vF.at(7), vF.at(3) } };
161  C.beTProductOf(F, F);
162  CC.beProductOf(C, C);
163  B.beProductTOf(F, F);
164  invF.beInverseOf(F);
165  FC.beProductOf(F, C);
166 
167  J = F.giveDeterminant();
168  lnJ = log(J);
169  I1 = C.at(1, 1) + C.at(2, 2) + C.at(3, 3);
170 
171  I2 = 0.5 * ( I1 * I1 - CC.at(1, 1) - CC.at(2, 2) - CC.at(3, 3) );
172 
173 
174  answer.at(1, 1) = ( 2. * C1 * ( 5. * I1 * invF.at(1, 1) * invF.at(1, 1) - 12. * F.at(1, 1) * invF.at(1, 1) + 9. ) ) / ( 9. * pow(J, 2. / 3.) );
175  answer.at(1, 2) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
176  answer.at(1, 3) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
177  answer.at(1, 4) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 2) - 3. * I1 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
178  answer.at(1, 5) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 1) + 6. * F.at(1, 3) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
179  answer.at(1, 6) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 1) + 6. * F.at(1, 2) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
180  answer.at(1, 7) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 3) - 3. * I1 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
181  answer.at(1, 8) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
182  answer.at(1, 9) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
183 
184 
185 
186 
187  answer.at(2, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
188  answer.at(2, 2) = ( 2. * C1 * ( 5. * I1 * invF.at(2, 2) * invF.at(2, 2) - 12. * F.at(2, 2) * invF.at(2, 2) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
189  answer.at(2, 3) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 3) - 3. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
190  answer.at(2, 4) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(2, 2) - 5. * I1 * invF.at(2, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
191  answer.at(2, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(3, 1) - 3. * I1 * invF.at(2, 1) * invF.at(3, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
192  answer.at(2, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
193  answer.at(2, 7) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(2, 2) - 5. * I1 * invF.at(2, 2) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
194  answer.at(2, 8) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 3) - 2. * I1 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
195  answer.at(2, 9) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
196 
197 
198  answer.at(3, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
199  answer.at(3, 2) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 3) - 3. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
200  answer.at(3, 3) = ( 2. * C1 * ( 5. * I1 * invF.at(3, 3) * invF.at(3, 3) - 12. * F.at(3, 3) * invF.at(3, 3) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
201  answer.at(3, 4) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(3, 2) - 5. * I1 * invF.at(3, 2) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
202  answer.at(3, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(3, 1) - 5. * I1 * invF.at(3, 1) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
203  answer.at(3, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 3) - 3. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
204  answer.at(3, 7) = -( 2. * C1 * ( 6. * F.at(3, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 3) - 5. * I1 * invF.at(2, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
205  answer.at(3, 8) = -( 2. * C1 * ( 6. * F.at(3, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 3) - 5. * I1 * invF.at(1, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
206  answer.at(3, 9) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
207 
208 
209 
210 
211  answer.at(4, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 2) - 3. * I1 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
212  answer.at(4, 2) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(2, 2) - 5. * I1 * invF.at(2, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
213  answer.at(4, 3) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(3, 2) - 5. * I1 * invF.at(3, 2) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
214  answer.at(4, 4) = ( 2. * C1 * ( 5. * I1 * invF.at(3, 2) * invF.at(3, 2) - 12. * F.at(2, 3) * invF.at(3, 2) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
215  answer.at(4, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(3, 1) - 5. * I1 * invF.at(3, 1) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
216  answer.at(4, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 2) - 3. * I1 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
217  answer.at(4, 7) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(3, 2) - 3. * I1 * invF.at(2, 2) * invF.at(3, 3) - 2. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
218  answer.at(4, 8) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(3, 2) - 3. * I1 * invF.at(1, 2) * invF.at(3, 3) - 2. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
219  answer.at(4, 9) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
220 
221 
222 
223  answer.at(5, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 1) + 6. * F.at(1, 3) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
224  answer.at(5, 2) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(3, 1) - 3. * I1 * invF.at(2, 1) * invF.at(3, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
225  answer.at(5, 3) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(3, 1) - 5. * I1 * invF.at(3, 1) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
226  answer.at(5, 4) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(3, 1) - 5. * I1 * invF.at(3, 1) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
227  answer.at(5, 5) = ( 2. * C1 * ( 5. * I1 * invF.at(3, 1) * invF.at(3, 1) - 12. * F.at(1, 3) * invF.at(3, 1) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
228  answer.at(5, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 1) + 6. * F.at(1, 3) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
229  answer.at(5, 7) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(3, 1) - 3. * I1 * invF.at(2, 1) * invF.at(3, 3) - 2. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
230  answer.at(5, 8) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(3, 1) - 3. * I1 * invF.at(1, 1) * invF.at(3, 3) - 2. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
231  answer.at(5, 9) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(3, 1) - 3. * I1 * invF.at(1, 1) * invF.at(3, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
232 
233 
234 
235  answer.at(6, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 1) + 6. * F.at(1, 2) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
236  answer.at(6, 2) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
237  answer.at(6, 3) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 3) - 3. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
238  answer.at(6, 4) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 2) - 3. * I1 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
239  answer.at(6, 5) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 1) + 6. * F.at(1, 3) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
240  answer.at(6, 6) = ( 2. * C1 * ( 5. * I1 * invF.at(2, 1) * invF.at(2, 1) - 12. * F.at(1, 2) * invF.at(2, 1) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
241  answer.at(6, 7) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
242  answer.at(6, 8) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 3) - 2. * I1 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
243  answer.at(6, 9) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
244 
245 
246 
247 
248  answer.at(7, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 3) - 3. * I1 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
249  answer.at(7, 2) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(2, 2) - 5. * I1 * invF.at(2, 2) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
250  answer.at(7, 3) = -( 2. * C1 * ( 6. * F.at(3, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 3) - 5. * I1 * invF.at(2, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
251  answer.at(7, 4) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(3, 2) - 3. * I1 * invF.at(2, 2) * invF.at(3, 3) - 2. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
252  answer.at(7, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(3, 1) - 3. * I1 * invF.at(2, 1) * invF.at(3, 3) - 2. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
253  answer.at(7, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
254  answer.at(7, 7) = ( 2. * C1 * ( 5. * I1 * invF.at(2, 3) * invF.at(2, 3) - 12. * F.at(3, 2) * invF.at(2, 3) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
255  answer.at(7, 8) = -( 2. * C1 * ( 6. * F.at(3, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 3) - 5. * I1 * invF.at(1, 3) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
256  answer.at(7, 9) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 3) - 3. * I1 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
257 
258 
259 
260 
261 
262  answer.at(8, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
263  answer.at(8, 2) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 3) - 2. * I1 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
264  answer.at(8, 3) = -( 2. * C1 * ( 6. * F.at(3, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 3) - 5. * I1 * invF.at(1, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
265  answer.at(8, 4) = -( 2. * C1 * ( 6. * F.at(2, 3) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(3, 2) - 3. * I1 * invF.at(1, 2) * invF.at(3, 3) - 2. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
266  answer.at(8, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(3, 1) - 3. * I1 * invF.at(1, 1) * invF.at(3, 3) - 2. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
267  answer.at(8, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 3) - 2. * I1 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
268  answer.at(8, 7) = -( 2. * C1 * ( 6. * F.at(3, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 3) - 5. * I1 * invF.at(1, 3) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
269  answer.at(8, 8) = ( 2. * C1 * ( 5. * I1 * invF.at(1, 3) * invF.at(1, 3) - 12. * F.at(3, 1) * invF.at(1, 3) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
270  answer.at(8, 9) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
271 
272 
273 
274 
275 
276  answer.at(9, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
277  answer.at(9, 2) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
278  answer.at(9, 3) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
279  answer.at(9, 4) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 2) + 6. * F.at(2, 3) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
280  answer.at(9, 5) = -( 2. * C1 * ( 6. * F.at(1, 3) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(3, 1) - 3. * I1 * invF.at(1, 1) * invF.at(3, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
281  answer.at(9, 6) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
282  answer.at(9, 7) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 3) + 6. * F.at(3, 2) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 3) - 3. * I1 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
283  answer.at(9, 8) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(1, 3) + 6. * F.at(3, 1) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(1, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
284  answer.at(9, 9) = ( 2. * C1 * ( 5. * I1 * invF.at(1, 2) * invF.at(1, 2) - 12. * F.at(2, 1) * invF.at(1, 2) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
285 
286 
287 
289 
290 
291  answer.at(1, 1) = answer.at(1, 1) + ( 2. * C2 * ( 9. * F.at(1, 1) * F.at(1, 1) - 24. * I1 * F.at(1, 1) * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 1) + 24. * FC.at(1, 1) * invF.at(1, 1) + 9. * B.at(1, 1) - 9. * C.at(1, 1) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
292  answer.at(1, 2) = answer.at(1, 2) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 1) + 18. * F.at(1, 1) * F.at(2, 2) - 9. * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(1, 1) - 8. * I2 * invF.at(1, 1) * invF.at(2, 2) + 6. * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
293  answer.at(1, 3) = answer.at(1, 3) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 1) - 8. * I2 * invF.at(1, 1) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
294  answer.at(1, 4) = answer.at(1, 4) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 3) + 12. * FC.at(2, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(2, 3) - 9. * F.at(1, 3) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(1, 1) - 8. * I2 * invF.at(1, 1) * invF.at(3, 2) + 6 * I2 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
295  answer.at(1, 5) = answer.at(1, 5) - ( 2. * C2 * ( 9. * C.at(3, 1) - 12. * FC.at(1, 1) * invF.at(1, 3) - 12. * FC.at(1, 3) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 3) + 12. * F.at(1, 1) * I1 * invF.at(3, 1) + 12. * F.at(1, 3) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
296  answer.at(1, 6) = answer.at(1, 6) - ( 2. * C2 * ( 9. * C.at(2, 1) - 12. * FC.at(1, 1) * invF.at(1, 2) - 12. * FC.at(1, 2) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 2) + 12. * F.at(1, 1) * I1 * invF.at(2, 1) + 12. * F.at(1, 2) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
297  answer.at(1, 7) = answer.at(1, 7) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 2) - 9. * F.at(1, 2) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(2, 3) + 6 * I2 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
298  answer.at(1, 8) = answer.at(1, 8) + ( 2. * C2 * ( 9. * B.at(1, 3) + 12. * FC.at(1, 1) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
299  answer.at(1, 9) = answer.at(1, 9) + ( 2. * C2 * ( 9. * B.at(1, 2) + 12. * FC.at(1, 1) * invF.at(2, 1) + 12. * FC.at(2, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
300 
301 
302 
303 
304 
305  answer.at(2, 1) = answer.at(2, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(2, 2) - 9. * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(2, 2) + 6 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
306  answer.at(2, 2) = answer.at(2, 2) + ( 2. * C2 * ( 9. * F.at(2, 2) * F.at(2, 2) - 24. * I1 * F.at(2, 2) * invF.at(2, 2) - 2. * I2 * invF.at(2, 2) * invF.at(2, 2) + 24. * FC.at(2, 2) * invF.at(2, 2) + 9. * B.at(2, 2) - 9. * C.at(2, 2) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
307  answer.at(2, 3) = answer.at(2, 3) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 2) + 18 * F.at(2, 2) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 3) + 6 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
308  answer.at(2, 4) = answer.at(2, 4) - ( 2. * C2 * ( 9. * C.at(3, 2) - 12. * FC.at(2, 2) * invF.at(2, 3) - 12. * FC.at(2, 3) * invF.at(2, 2) - 9. * F.at(2, 2) * F.at(2, 3) + 12. * F.at(2, 2) * I1 * invF.at(3, 2) + 12. * F.at(2, 3) * I1 * invF.at(2, 2) + 2. * I2 * invF.at(2, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
309  answer.at(2, 5) = answer.at(2, 5) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 3) - 9. * F.at(1, 2) * F.at(2, 3) + 18 * F.at(1, 3) * F.at(2, 2) - 12. * F.at(1, 3) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(2, 1) * invF.at(3, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
310  answer.at(2, 6) = answer.at(2, 6) + ( 2. * C2 * ( 9. * B.at(2, 1) + 12. * FC.at(1, 2) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 2) + 9. * F.at(1, 2) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
311  answer.at(2, 7) = answer.at(2, 7) + ( 2. * C2 * ( 9. * B.at(2, 3) + 12. * FC.at(2, 2) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(2, 2) + 9. * F.at(2, 2) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(2, 2) - 2. * I2 * invF.at(2, 2) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
312  answer.at(2, 8) = answer.at(2, 8) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(2, 2) - 9. * F.at(2, 1) * F.at(3, 2) + 18 * F.at(2, 2) * F.at(3, 1) - 12. * F.at(2, 2) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(2, 2) + 6 * I2 * invF.at(1, 2) * invF.at(2, 3) - 8 * I2 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
313  answer.at(2, 9) = answer.at(2, 9) - ( 2. * C2 * ( 9. * C.at(1, 2) - 12. * FC.at(2, 1) * invF.at(2, 2) - 12. * FC.at(2, 2) * invF.at(2, 1) - 9. * F.at(2, 1) * F.at(2, 2) + 12. * F.at(2, 1) * I1 * invF.at(2, 2) + 12. * F.at(2, 2) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
314 
315 
316 
317 
318 
319  answer.at(3, 1) = answer.at(3, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
320  answer.at(3, 2) = answer.at(3, 2) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 2) + 18 * F.at(2, 2) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 3) + 6 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
321  answer.at(3, 3) = answer.at(3, 3) + ( 2. * C2 * ( 9. * F.at(3, 3) * F.at(3, 3) - 24. * I1 * F.at(3, 3) * invF.at(3, 3) - 2. * I2 * invF.at(3, 3) * invF.at(3, 3) + 24. * FC.at(3, 3) * invF.at(3, 3) + 9. * B.at(3, 3) - 9. * C.at(3, 3) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
322  answer.at(3, 4) = answer.at(3, 4) + ( 2. * C2 * ( 9. * B.at(3, 2) + 12. * FC.at(2, 3) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 3) + 9. * F.at(2, 3) * F.at(3, 3) - 12. * F.at(2, 3) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(3, 2) - 2. * I2 * invF.at(3, 2) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
323  answer.at(3, 5) = answer.at(3, 5) + ( 2. * C2 * ( 9. * B.at(3, 1) + 12. * FC.at(1, 3) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 3) + 9. * F.at(1, 3) * F.at(3, 3) - 12. * F.at(1, 3) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(3, 1) - 2. * I2 * invF.at(3, 1) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
324  answer.at(3, 6) = answer.at(3, 6) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 2) + 18 * F.at(1, 2) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 3) + 6 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
325  answer.at(3, 7) = answer.at(3, 7) - ( 2. * C2 * ( 9. * C.at(2, 3) - 12. * FC.at(3, 2) * invF.at(3, 3) - 12. * FC.at(3, 3) * invF.at(3, 2) - 9. * F.at(3, 2) * F.at(3, 3) + 12. * F.at(3, 2) * I1 * invF.at(3, 3) + 12. * F.at(3, 3) * I1 * invF.at(2, 3) + 2. * I2 * invF.at(2, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
326  answer.at(3, 8) = answer.at(3, 8) - ( 2. * C2 * ( 9. * C.at(1, 3) - 12. * FC.at(3, 1) * invF.at(3, 3) - 12. * FC.at(3, 3) * invF.at(3, 1) - 9. * F.at(3, 1) * F.at(3, 3) + 12. * F.at(3, 1) * I1 * invF.at(3, 3) + 12. * F.at(3, 3) * I1 * invF.at(1, 3) + 2. * I2 * invF.at(1, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
327  answer.at(3, 9) = answer.at(3, 9) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 1) + 18 * F.at(2, 1) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
328 
329 
330 
331  answer.at(4, 1) = answer.at(4, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 2) + 12. * FC.at(2, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(2, 3) - 9. * F.at(1, 3) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(3, 2) + 6 * I2 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
332  answer.at(4, 2) = answer.at(4, 2) - ( 2. * C2 * ( 9. * C.at(2, 3) - 12. * FC.at(2, 2) * invF.at(3, 2) - 12. * FC.at(2, 3) * invF.at(2, 2) - 9. * F.at(2, 2) * F.at(2, 3) + 12. * F.at(2, 2) * I1 * invF.at(3, 2) + 12. * F.at(2, 3) * I1 * invF.at(2, 2) + 2. * I2 * invF.at(2, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
333  answer.at(4, 3) = answer.at(4, 3) + ( 2. * C2 * ( 9. * B.at(2, 3) + 12. * FC.at(2, 3) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(3, 2) + 9. * F.at(2, 3) * F.at(3, 3) - 12. * F.at(2, 3) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(3, 2) - 2. * I2 * invF.at(3, 2) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
334  answer.at(4, 4) = answer.at(4, 4) + ( 2. * C2 * ( 9. * F.at(2, 3) * F.at(2, 3) - 24. * I1 * F.at(2, 3) * invF.at(3, 2) - 2. * I2 * invF.at(3, 2) * invF.at(3, 2) + 12. * FC.at(2, 3) * invF.at(3, 2) + 9. * B.at(2, 2) - 9. * C.at(3, 3) + 9. * I1 + 12. * FC.at(2, 3) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
335  answer.at(4, 5) = answer.at(4, 5) + ( 2. * C2 * ( 9. * B.at(2, 1) + 12. * FC.at(1, 3) * invF.at(3, 2) + 12. * FC.at(2, 3) * invF.at(1, 3) + 9. * F.at(1, 3) * F.at(2, 3) - 12. * F.at(1, 3) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(3, 1) - 2. * I2 * invF.at(3, 1) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
336  answer.at(4, 6) = answer.at(4, 6) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 2) + 12. * FC.at(2, 3) * invF.at(1, 2) + 18 * F.at(1, 2) * F.at(2, 3) - 9. * F.at(1, 3) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 2) + 6 * I2 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
337  answer.at(4, 7) = answer.at(4, 7) + ( 2. * C2 * ( 12. * FC.at(2, 3) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(3, 2) - 9. * F.at(2, 2) * F.at(3, 3) + 18 * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(3, 2) + 6 * I2 * invF.at(2, 2) * invF.at(3, 3) - 8 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
338  answer.at(4, 8) = answer.at(4, 8) + ( 2. * C2 * ( 12. * FC.at(2, 3) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(3, 2) - 9. * F.at(2, 1) * F.at(3, 3) + 18 * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 3) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(3, 2) + 6 * I2 * invF.at(1, 2) * invF.at(3, 3) - 8 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
339  answer.at(4, 9) = answer.at(4, 9) - ( 2. * C2 * ( 9. * C.at(1, 3) - 12. * FC.at(2, 1) * invF.at(3, 2) - 12. * FC.at(2, 3) * invF.at(2, 1) - 9. * F.at(2, 1) * F.at(2, 3) + 12. * F.at(2, 1) * I1 * invF.at(3, 2) + 12. * F.at(2, 3) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
340 
341 
342  answer.at(5, 1) = answer.at(5, 1) - ( 2. * C2 * ( 9. * C.at(1, 3) - 12. * FC.at(1, 1) * invF.at(3, 1) - 12. * FC.at(1, 3) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 3) + 12. * F.at(1, 1) * I1 * invF.at(3, 1) + 12. * F.at(1, 3) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
343  answer.at(5, 2) = answer.at(5, 2) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(3, 1) - 9. * F.at(1, 2) * F.at(2, 3) + 18 * F.at(1, 3) * F.at(2, 2) - 12. * F.at(1, 3) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(2, 1) * invF.at(3, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
344  answer.at(5, 3) = answer.at(5, 3) + ( 2. * C2 * ( 9. * B.at(1, 3) + 12. * FC.at(1, 3) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(3, 1) + 9. * F.at(1, 3) * F.at(3, 3) - 12. * F.at(1, 3) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(3, 1) - 2. * I2 * invF.at(3, 1) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
345  answer.at(5, 4) = answer.at(5, 4) + ( 2. * C2 * ( 9. * B.at(1, 2) + 12. * FC.at(1, 3) * invF.at(2, 3) + 12. * FC.at(2, 3) * invF.at(3, 1) + 9. * F.at(1, 3) * F.at(2, 3) - 12. * F.at(1, 3) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(3, 1) - 2. * I2 * invF.at(3, 1) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
346  answer.at(5, 5) = answer.at(5, 5) + ( 2. * C2 * ( 9. * F.at(1, 3) * F.at(1, 3) - 24. * I1 * F.at(1, 3) * invF.at(3, 1) - 2. * I2 * invF.at(3, 1) * invF.at(3, 1) + 12. * FC.at(1, 3) * invF.at(3, 1) + 9. * B.at(1, 1) - 9. * C.at(3, 3) + 9. * I1 + 12. * FC.at(1, 3) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
347  answer.at(5, 6) = answer.at(5, 6) - ( 2. * C2 * ( 9. * C.at(2, 3) - 12. * FC.at(1, 2) * invF.at(3, 1) - 12. * FC.at(1, 3) * invF.at(1, 2) - 9. * F.at(1, 2) * F.at(1, 3) + 12. * F.at(1, 2) * I1 * invF.at(3, 1) + 12. * F.at(1, 3) * I1 * invF.at(2, 1) + 2. * I2 * invF.at(2, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
348  answer.at(5, 7) = answer.at(5, 7) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(3, 1) - 9. * F.at(1, 2) * F.at(3, 3) + 18 * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 3) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(2, 1) * invF.at(3, 3) - 8 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
349  answer.at(5, 8) = answer.at(5, 8) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(3, 1) - 9. * F.at(1, 1) * F.at(3, 3) + 18 * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 3) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(1, 1) * invF.at(3, 3) - 8 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
350  answer.at(5, 9) = answer.at(5, 9) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(2, 1) + 12. * FC.at(2, 1) * invF.at(3, 1) - 9. * F.at(1, 1) * F.at(2, 3) + 18 * F.at(1, 3) * F.at(2, 1) - 12. * F.at(1, 3) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(1, 1) * invF.at(3, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
351 
352 
353 
354  answer.at(6, 1) = answer.at(6, 1) - ( 2. * C2 * ( 9. * C.at(1, 2) - 12. * FC.at(1, 1) * invF.at(2, 1) - 12. * FC.at(1, 2) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 2) + 12. * F.at(1, 1) * I1 * invF.at(2, 1) + 12. * F.at(1, 2) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
355  answer.at(6, 2) = answer.at(6, 2) + ( 2. * C2 * ( 9. * B.at(1, 2) + 12. * FC.at(1, 2) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(2, 1) + 9. * F.at(1, 2) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
356  answer.at(6, 3) = answer.at(6, 3) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 1) + 18 * F.at(1, 2) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 3) + 6 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
357  answer.at(6, 4) = answer.at(6, 4) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(2, 3) + 12. * FC.at(2, 3) * invF.at(2, 1) + 18 * F.at(1, 2) * F.at(2, 3) - 9. * F.at(1, 3) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 2) + 6 * I2 * invF.at(2, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
358  answer.at(6, 5) = answer.at(6, 5) - ( 2. * C2 * ( 9. * C.at(3, 2) - 12. * FC.at(1, 2) * invF.at(1, 3) - 12. * FC.at(1, 3) * invF.at(2, 1) - 9. * F.at(1, 2) * F.at(1, 3) + 12. * F.at(1, 2) * I1 * invF.at(3, 1) + 12. * F.at(1, 3) * I1 * invF.at(2, 1) + 2. * I2 * invF.at(2, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
359  answer.at(6, 6) = answer.at(6, 6) + ( 2. * C2 * ( 9. * F.at(1, 2) * F.at(1, 2) - 24. * I1 * F.at(1, 2) * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 1) + 12. * FC.at(1, 2) * invF.at(2, 1) + 9. * B.at(1, 1) - 9. * C.at(2, 2) + 9. * I1 + 12. * FC.at(1, 2) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
360  answer.at(6, 7) = answer.at(6, 7) + ( 2. * C2 * ( 9. * B.at(1, 3) + 12. * FC.at(1, 2) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(2, 1) + 9. * F.at(1, 2) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
361  answer.at(6, 8) = answer.at(6, 8) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(2, 1) - 9. * F.at(1, 1) * F.at(3, 2) + 18 * F.at(1, 2) * F.at(3, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 3) - 8 * I2 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
362  answer.at(6, 9) = answer.at(6, 9) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(2, 1) + 12. * FC.at(2, 1) * invF.at(2, 1) - 9. * F.at(1, 1) * F.at(2, 2) + 18 * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
363 
364 
365 
366 
367  answer.at(7, 1) = answer.at(7, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 2) - 9. * F.at(1, 2) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(2, 3) + 6 * I2 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
368  answer.at(7, 2) = answer.at(7, 2) + ( 2. * C2 * ( 9. * B.at(3, 2) + 12. * FC.at(2, 2) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(2, 2) + 9. * F.at(2, 2) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(2, 2) - 2. * I2 * invF.at(2, 2) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
369  answer.at(7, 3) = answer.at(7, 3) - ( 2. * C2 * ( 9. * C.at(3, 2) - 12. * FC.at(3, 2) * invF.at(3, 3) - 12. * FC.at(3, 3) * invF.at(2, 3) - 9. * F.at(3, 2) * F.at(3, 3) + 12. * F.at(3, 2) * I1 * invF.at(3, 3) + 12. * F.at(3, 3) * I1 * invF.at(2, 3) + 2. * I2 * invF.at(2, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
370  answer.at(7, 4) = answer.at(7, 4) + ( 2. * C2 * ( 12. * FC.at(2, 3) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(2, 3) - 9. * F.at(2, 2) * F.at(3, 3) + 18 * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 3) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(3, 2) + 6 * I2 * invF.at(2, 2) * invF.at(3, 3) - 8 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
371  answer.at(7, 5) = answer.at(7, 5) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(1, 3) - 9. * F.at(1, 2) * F.at(3, 3) + 18 * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 3) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(2, 1) * invF.at(3, 3) - 8 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
372  answer.at(7, 6) = answer.at(7, 6) + ( 2. * C2 * ( 9. * B.at(3, 1) + 12. * FC.at(1, 2) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(1, 2) + 9. * F.at(1, 2) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
373  answer.at(7, 7) = answer.at(7, 7) + ( 2. * C2 * ( 9. * F.at(3, 2) * F.at(3, 2) - 24. * I1 * F.at(3, 2) * invF.at(2, 3) - 2. * I2 * invF.at(2, 3) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(2, 3) + 9. * B.at(3, 3) - 9. * C.at(2, 2) + 9. * I1 + 12. * FC.at(3, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
374  answer.at(7, 8) = answer.at(7, 8) - ( 2. * C2 * ( 9. * C.at(1, 2) - 12. * FC.at(3, 1) * invF.at(2, 3) - 12. * FC.at(3, 2) * invF.at(3, 1) - 9. * F.at(3, 1) * F.at(3, 2) + 12. * F.at(3, 1) * I1 * invF.at(2, 3) + 12. * F.at(3, 2) * I1 * invF.at(1, 3) + 2. * I2 * invF.at(1, 3) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
375  answer.at(7, 9) = answer.at(7, 9) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(2, 3) + 12. * FC.at(3, 2) * invF.at(2, 1) + 18 * F.at(2, 1) * F.at(3, 2) - 9. * F.at(2, 2) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 3) + 6 * I2 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
376 
377 
378 
379  answer.at(8, 1) = answer.at(8, 1) + ( 2. * C2 * ( 9. * B.at(3, 1) + 12. * FC.at(1, 1) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
380  answer.at(8, 2) = answer.at(8, 2) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(2, 2) - 9. * F.at(2, 1) * F.at(3, 2) + 18 * F.at(2, 2) * F.at(3, 1) - 12. * F.at(2, 2) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(2, 2) + 6 * I2 * invF.at(1, 2) * invF.at(2, 3) - 8 * I2 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
381  answer.at(8, 3) = answer.at(8, 3) - ( 2. * C2 * ( 9. * C.at(3, 1) - 12. * FC.at(3, 1) * invF.at(3, 3) - 12. * FC.at(3, 3) * invF.at(1, 3) - 9. * F.at(3, 1) * F.at(3, 3) + 12. * F.at(3, 1) * I1 * invF.at(3, 3) + 12. * F.at(3, 3) * I1 * invF.at(1, 3) + 2. * I2 * invF.at(1, 3) * invF.at(3, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
382  answer.at(8, 4) = answer.at(8, 4) + ( 2. * C2 * ( 12. * FC.at(2, 3) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(2, 3) - 9. * F.at(2, 1) * F.at(3, 3) + 18 * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 3) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(3, 2) + 6 * I2 * invF.at(1, 2) * invF.at(3, 3) - 8 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
383  answer.at(8, 5) = answer.at(8, 5) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(1, 3) - 9. * F.at(1, 1) * F.at(3, 3) + 18 * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 3) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(1, 1) * invF.at(3, 3) - 8 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
384  answer.at(8, 6) = answer.at(8, 6) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(1, 2) - 9. * F.at(1, 1) * F.at(3, 2) + 18 * F.at(1, 2) * F.at(3, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 3) - 8 * I2 * invF.at(2, 1) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
385  answer.at(8, 7) = answer.at(8, 7) - ( 2. * C2 * ( 9. * C.at(2, 1) - 12. * FC.at(3, 1) * invF.at(3, 2) - 12. * FC.at(3, 2) * invF.at(1, 3) - 9. * F.at(3, 1) * F.at(3, 2) + 12. * F.at(3, 1) * I1 * invF.at(2, 3) + 12. * F.at(3, 2) * I1 * invF.at(1, 3) + 2. * I2 * invF.at(1, 3) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
386  answer.at(8, 8) = answer.at(8, 8) + ( 2. * C2 * ( 9. * F.at(3, 1) * F.at(3, 1) - 24. * I1 * F.at(3, 1) * invF.at(1, 3) - 2. * I2 * invF.at(1, 3) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(1, 3) + 9. * B.at(3, 3) - 9. * C.at(1, 1) + 9. * I1 + 12. * FC.at(3, 1) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
387  answer.at(8, 9) = answer.at(8, 9) + ( 2. * C2 * ( 9. * B.at(3, 2) + 12. * FC.at(2, 1) * invF.at(1, 3) + 12. * FC.at(3, 1) * invF.at(2, 1) + 9. * F.at(2, 1) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(1, 2) - 2. * I2 * invF.at(1, 2) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
388 
389 
390 
391 
392  answer.at(9, 1) = answer.at(9, 1) + ( 2. * C2 * ( 9. * B.at(2, 1) + 12. * FC.at(1, 1) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
393  answer.at(9, 2) = answer.at(9, 2) - ( 2. * C2 * ( 9. * C.at(2, 1) - 12. * FC.at(2, 1) * invF.at(2, 2) - 12. * FC.at(2, 2) * invF.at(1, 2) - 9. * F.at(2, 1) * F.at(2, 2) + 12. * F.at(2, 1) * I1 * invF.at(2, 2) + 12. * F.at(2, 2) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
394  answer.at(9, 3) = answer.at(9, 3) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 2) + 18 * F.at(2, 1) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
395  answer.at(9, 4) = answer.at(9, 4) - ( 2. * C2 * ( 9. * C.at(3, 1) - 12. * FC.at(2, 1) * invF.at(2, 3) - 12. * FC.at(2, 3) * invF.at(1, 2) - 9. * F.at(2, 1) * F.at(2, 3) + 12. * F.at(2, 1) * I1 * invF.at(3, 2) + 12. * F.at(2, 3) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
396  answer.at(9, 5) = answer.at(9, 5) + ( 2. * C2 * ( 12. * FC.at(1, 3) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 3) - 9. * F.at(1, 1) * F.at(2, 3) + 18 * F.at(1, 3) * F.at(2, 1) - 12. * F.at(1, 3) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(3, 1) + 6 * I2 * invF.at(1, 1) * invF.at(3, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
397  answer.at(9, 6) = answer.at(9, 6) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 2) - 9. * F.at(1, 1) * F.at(2, 2) + 18 * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
398  answer.at(9, 7) = answer.at(9, 7) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(3, 2) + 12. * FC.at(3, 2) * invF.at(1, 2) + 18 * F.at(2, 1) * F.at(3, 2) - 9. * F.at(2, 2) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(2, 3) - 12. * F.at(3, 2) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 3) + 6 * I2 * invF.at(1, 3) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
399  answer.at(9, 8) = answer.at(9, 8) + ( 2. * C2 * ( 9. * B.at(2, 3) + 12. * FC.at(2, 1) * invF.at(3, 1) + 12. * FC.at(3, 1) * invF.at(1, 2) + 9. * F.at(2, 1) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(1, 3) - 12. * F.at(3, 1) * I1 * invF.at(1, 2) - 2. * I2 * invF.at(1, 2) * invF.at(1, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
400  answer.at(9, 9) = answer.at(9, 9) + ( 2. * C2 * ( 9. * F.at(2, 1) * F.at(2, 1) - 24. * I1 * F.at(2, 1) * invF.at(1, 2) - 2. * I2 * invF.at(1, 2) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 2) + 9. * B.at(2, 2) - 9. * C.at(1, 1) + 9. * I1 + 12. * FC.at(2, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
401 
403 
404 
405 
406 
407 
408  answer.at(1, 1) = answer.at(1, 1) - K *invF.at(1, 1) * invF.at(1, 1) * ( lnJ - 1. );
409  answer.at(1, 2) = answer.at(1, 2) + K * ( invF.at(1, 1) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 1) * lnJ );
410  answer.at(1, 3) = answer.at(1, 3) + K * ( invF.at(1, 1) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 1) * lnJ );
411  answer.at(1, 4) = answer.at(1, 4) + K * ( invF.at(1, 1) * invF.at(3, 2) - invF.at(1, 2) * invF.at(3, 1) * lnJ );
412  answer.at(1, 5) = answer.at(1, 5) - K *invF.at(1, 1) * invF.at(3, 1) * ( lnJ - 1. );
413  answer.at(1, 6) = answer.at(1, 6) - K *invF.at(1, 1) * invF.at(2, 1) * ( lnJ - 1. );
414  answer.at(1, 7) = answer.at(1, 7) + K * ( invF.at(1, 1) * invF.at(2, 3) - invF.at(2, 1) * invF.at(1, 3) * lnJ );
415  answer.at(1, 8) = answer.at(1, 8) - K *invF.at(1, 1) * invF.at(1, 3) * ( lnJ - 1. );
416  answer.at(1, 9) = answer.at(1, 9) - K *invF.at(1, 1) * invF.at(1, 2) * ( lnJ - 1. );
417 
418 
419 
420 
421  answer.at(2, 1) = answer.at(2, 1) + K * ( invF.at(1, 1) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 1) * lnJ );
422  answer.at(2, 2) = answer.at(2, 2) - K *invF.at(2, 2) * invF.at(2, 2) * ( lnJ - 1. );
423  answer.at(2, 3) = answer.at(2, 3) + K * ( invF.at(2, 2) * invF.at(3, 3) - invF.at(2, 3) * invF.at(3, 2) * lnJ );
424  answer.at(2, 4) = answer.at(2, 4) - K *invF.at(2, 2) * invF.at(3, 2) * ( lnJ - 1. );
425  answer.at(2, 5) = answer.at(2, 5) + K * ( invF.at(2, 2) * invF.at(3, 1) - invF.at(2, 1) * invF.at(3, 2) * lnJ );
426  answer.at(2, 6) = answer.at(2, 6) - K *invF.at(2, 1) * invF.at(2, 2) * ( lnJ - 1. );
427  answer.at(2, 7) = answer.at(2, 7) - K *invF.at(2, 2) * invF.at(2, 3) * ( lnJ - 1. );
428  answer.at(2, 8) = answer.at(2, 8) + K * ( invF.at(1, 3) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 3) * lnJ );
429  answer.at(2, 9) = answer.at(2, 9) - K *invF.at(1, 2) * invF.at(2, 2) * ( lnJ - 1. );
430 
431 
432 
433 
434  answer.at(3, 1) = answer.at(3, 1) + K * ( invF.at(1, 1) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 1) * lnJ );
435  answer.at(3, 2) = answer.at(3, 2) + K * ( invF.at(2, 2) * invF.at(3, 3) - invF.at(2, 3) * invF.at(3, 2) * lnJ );
436  answer.at(3, 3) = answer.at(3, 3) - K *invF.at(3, 3) * invF.at(3, 3) * ( lnJ - 1. );
437  answer.at(3, 4) = answer.at(3, 4) - K *invF.at(3, 2) * invF.at(3, 3) * ( lnJ - 1. );
438  answer.at(3, 5) = answer.at(3, 5) - K *invF.at(3, 1) * invF.at(3, 3) * ( lnJ - 1. );
439  answer.at(3, 6) = answer.at(3, 6) + K * ( invF.at(2, 1) * invF.at(3, 3) - invF.at(3, 1) * invF.at(2, 3) * lnJ );
440  answer.at(3, 7) = answer.at(3, 7) - K *invF.at(2, 3) * invF.at(3, 3) * ( lnJ - 1. );
441  answer.at(3, 8) = answer.at(3, 8) - K *invF.at(1, 3) * invF.at(3, 3) * ( lnJ - 1. );
442  answer.at(3, 9) = answer.at(3, 9) + K * ( invF.at(1, 2) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 2) * lnJ );
443 
444 
445 
446 
447  answer.at(4, 1) = answer.at(4, 1) + K * ( invF.at(1, 1) * invF.at(3, 2) - invF.at(1, 2) * invF.at(3, 1) * lnJ );
448  answer.at(4, 2) = answer.at(4, 2) - K *invF.at(2, 2) * invF.at(3, 2) * ( lnJ - 1. );
449  answer.at(4, 3) = answer.at(4, 3) - K *invF.at(3, 2) * invF.at(3, 3) * ( lnJ - 1. );
450  answer.at(4, 4) = answer.at(4, 4) - K *invF.at(3, 2) * invF.at(3, 2) * ( lnJ - 1. );
451  answer.at(4, 5) = answer.at(4, 5) - K *invF.at(3, 1) * invF.at(3, 2) * ( lnJ - 1. );
452  answer.at(4, 6) = answer.at(4, 6) + K * ( invF.at(2, 1) * invF.at(3, 2) - invF.at(2, 2) * invF.at(3, 1) * lnJ );
453  answer.at(4, 7) = answer.at(4, 7) + K * ( invF.at(2, 3) * invF.at(3, 2) - invF.at(2, 2) * invF.at(3, 3) * lnJ );
454  answer.at(4, 8) = answer.at(4, 8) + K * ( invF.at(1, 3) * invF.at(3, 2) - invF.at(1, 2) * invF.at(3, 3) * lnJ );
455  answer.at(4, 9) = answer.at(4, 9) - K *invF.at(1, 2) * invF.at(3, 2) * ( lnJ - 1. );
456 
457 
458 
459  answer.at(5, 1) = answer.at(5, 1) - K *invF.at(1, 1) * invF.at(3, 1) * ( lnJ - 1. );
460  answer.at(5, 2) = answer.at(5, 2) + K * ( invF.at(2, 2) * invF.at(3, 1) - invF.at(2, 1) * invF.at(3, 2) * lnJ );
461  answer.at(5, 3) = answer.at(5, 3) - K *invF.at(3, 1) * invF.at(3, 3) * ( lnJ - 1. );
462  answer.at(5, 4) = answer.at(5, 4) - K *invF.at(3, 1) * invF.at(3, 2) * ( lnJ - 1. );
463  answer.at(5, 5) = answer.at(5, 5) - K *invF.at(3, 1) * invF.at(3, 1) * ( lnJ - 1. );
464  answer.at(5, 6) = answer.at(5, 6) - K *invF.at(2, 1) * invF.at(3, 1) * ( lnJ - 1. );
465  answer.at(5, 7) = answer.at(5, 7) + K * ( invF.at(3, 1) * invF.at(2, 3) - invF.at(2, 1) * invF.at(3, 3) * lnJ );
466  answer.at(5, 8) = answer.at(5, 8) + K * ( invF.at(1, 3) * invF.at(3, 1) - invF.at(1, 1) * invF.at(3, 3) * lnJ );
467  answer.at(5, 9) = answer.at(5, 9) + K * ( invF.at(1, 2) * invF.at(3, 1) - invF.at(1, 1) * invF.at(3, 2) * lnJ );
468 
469 
470 
471 
472  answer.at(6, 1) = answer.at(6, 1) - K *invF.at(1, 1) * invF.at(2, 1) * ( lnJ - 1. );
473  answer.at(6, 2) = answer.at(6, 2) - K *invF.at(2, 1) * invF.at(2, 2) * ( lnJ - 1. );
474  answer.at(6, 3) = answer.at(6, 3) + K * ( invF.at(2, 1) * invF.at(3, 3) - invF.at(3, 1) * invF.at(2, 3) * lnJ );
475  answer.at(6, 4) = answer.at(6, 4) + K * ( invF.at(2, 1) * invF.at(3, 2) - invF.at(2, 2) * invF.at(3, 1) * lnJ );
476  answer.at(6, 5) = answer.at(6, 5) - K *invF.at(2, 1) * invF.at(3, 1) * ( lnJ - 1. );
477  answer.at(6, 6) = answer.at(6, 6) - K *invF.at(2, 1) * invF.at(2, 1) * ( lnJ - 1. );
478  answer.at(6, 7) = answer.at(6, 7) - K *invF.at(2, 1) * invF.at(2, 3) * ( lnJ - 1. );
479  answer.at(6, 8) = answer.at(6, 8) + K * ( invF.at(2, 1) * invF.at(1, 3) - invF.at(1, 1) * invF.at(2, 3) * lnJ );
480  answer.at(6, 9) = answer.at(6, 9) + K * ( invF.at(1, 2) * invF.at(2, 1) - invF.at(1, 1) * invF.at(2, 2) * lnJ );
481 
482 
483 
484 
485 
486  answer.at(7, 1) = answer.at(7, 1) + K * ( invF.at(1, 1) * invF.at(2, 3) - invF.at(2, 1) * invF.at(1, 3) * lnJ );
487  answer.at(7, 2) = answer.at(7, 2) - K *invF.at(2, 2) * invF.at(2, 3) * ( lnJ - 1. );
488  answer.at(7, 3) = answer.at(7, 3) - K *invF.at(2, 3) * invF.at(3, 3) * ( lnJ - 1. );
489  answer.at(7, 4) = answer.at(7, 4) + K * ( invF.at(2, 3) * invF.at(3, 2) - invF.at(2, 2) * invF.at(3, 3) * lnJ );
490  answer.at(7, 5) = answer.at(7, 5) + K * ( invF.at(3, 1) * invF.at(2, 3) - invF.at(2, 1) * invF.at(3, 3) * lnJ );
491  answer.at(7, 6) = answer.at(7, 6) - K *invF.at(2, 1) * invF.at(2, 3) * ( lnJ - 1. );
492  answer.at(7, 7) = answer.at(7, 7) - K *invF.at(2, 3) * invF.at(2, 3) * ( lnJ - 1. );
493  answer.at(7, 8) = answer.at(7, 8) - K *invF.at(1, 3) * invF.at(2, 3) * ( lnJ - 1. );
494  answer.at(7, 9) = answer.at(7, 9) + K * ( invF.at(1, 2) * invF.at(2, 3) - invF.at(1, 3) * invF.at(2, 2) * lnJ );
495 
496 
497 
498 
499  answer.at(8, 1) = answer.at(8, 1) - K *invF.at(1, 1) * invF.at(1, 3) * ( lnJ - 1. );
500  answer.at(8, 2) = answer.at(8, 2) + K * ( invF.at(1, 3) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 3) * lnJ );
501  answer.at(8, 3) = answer.at(8, 3) - K *invF.at(1, 3) * invF.at(3, 3) * ( lnJ - 1. );
502  answer.at(8, 4) = answer.at(8, 4) + K * ( invF.at(1, 3) * invF.at(3, 2) - invF.at(1, 2) * invF.at(3, 3) * lnJ );
503  answer.at(8, 5) = answer.at(8, 5) + K * ( invF.at(1, 3) * invF.at(3, 1) - invF.at(1, 1) * invF.at(3, 3) * lnJ );
504  answer.at(8, 6) = answer.at(8, 6) + K * ( invF.at(2, 1) * invF.at(1, 3) - invF.at(1, 1) * invF.at(2, 3) * lnJ );
505  answer.at(8, 7) = answer.at(8, 7) - K *invF.at(1, 3) * invF.at(2, 3) * ( lnJ - 1. );
506  answer.at(8, 8) = answer.at(8, 8) - K *invF.at(1, 3) * invF.at(1, 3) * ( lnJ - 1. );
507  answer.at(8, 9) = answer.at(8, 9) - K *invF.at(1, 2) * invF.at(1, 3) * ( lnJ - 1. );
508 
509 
510 
511 
512 
513  answer.at(9, 1) = answer.at(9, 1) - K *invF.at(1, 1) * invF.at(1, 2) * ( lnJ - 1. );
514  answer.at(9, 2) = answer.at(9, 2) - K *invF.at(1, 2) * invF.at(2, 2) * ( lnJ - 1. );
515  answer.at(9, 3) = answer.at(9, 3) + K * ( invF.at(1, 2) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 2) * lnJ );
516  answer.at(9, 4) = answer.at(9, 4) - K *invF.at(1, 2) * invF.at(3, 2) * ( lnJ - 1. );
517  answer.at(9, 5) = answer.at(9, 5) + K * ( invF.at(1, 2) * invF.at(3, 1) - invF.at(1, 1) * invF.at(3, 2) * lnJ );
518  answer.at(9, 6) = answer.at(9, 6) + K * ( invF.at(1, 2) * invF.at(2, 1) - invF.at(1, 1) * invF.at(2, 2) * lnJ );
519  answer.at(9, 7) = answer.at(9, 7) + K * ( invF.at(1, 2) * invF.at(2, 3) - invF.at(1, 3) * invF.at(2, 2) * lnJ );
520  answer.at(9, 8) = answer.at(9, 8) - K *invF.at(1, 2) * invF.at(1, 3) * ( lnJ - 1. );
521  answer.at(9, 9) = answer.at(9, 9) - K *invF.at(1, 2) * invF.at(1, 2) * ( lnJ - 1. );
522 }
523 
524 
525 
526 void
527 MooneyRivlinMaterial :: givePlaneStrainStiffMtrx_dPdF(FloatMatrix &answer,
528  MatResponseMode mode,
529  GaussPoint *gp, TimeStep *tStep)
530 {
531  StructuralMaterialStatus *status = static_cast< StructuralMaterialStatus * >( this->giveStatus(gp) );
532 
533  double I1, I2, J, lnJ;
534  FloatMatrix C, CC, invF, FC, B;
535 
536 
537  FloatArray vF = status->giveTempFVector();
538  FloatMatrix F = { { vF.at(1), vF.at(5), 0 }, { vF.at(4), vF.at(2), 0 }, { 0, 0, vF.at(3) } };
539  C.beTProductOf(F, F);
540  CC.beProductOf(C, C);
541  B.beProductTOf(F, F);
542  invF.beInverseOf(F);
543  FC.beProductOf(F, C);
544 
545  J = F.giveDeterminant();
546  lnJ = log(J);
547  I1 = C.at(1, 1) + C.at(2, 2) + C.at(3, 3);
548 
549  I2 = 0.5 * ( I1 * I1 - CC.at(1, 1) - CC.at(2, 2) - CC.at(3, 3) );
550 
551  answer.resize(5, 5);
552 
553  answer.at(1, 1) = ( 2. * C1 * ( 5. * I1 * invF.at(1, 1) * invF.at(1, 1) - 12. * F.at(1, 1) * invF.at(1, 1) + 9. ) ) / ( 9. * pow(J, 2. / 3.) );
554  answer.at(1, 2) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
555  answer.at(1, 3) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
556  answer.at(1, 4) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 1) + 6. * F.at(1, 2) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
557  answer.at(1, 5) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
558 
559  answer.at(2, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(2, 2) - 3. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
560  answer.at(2, 2) = ( 2. * C1 * ( 5. * I1 * invF.at(2, 2) * invF.at(2, 2) - 12. * F.at(2, 2) * invF.at(2, 2) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
561  answer.at(2, 3) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 3) - 3. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
562  answer.at(2, 4) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
563  answer.at(2, 5) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
564 
565 
566  answer.at(3, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 1) - 2. * I1 * invF.at(1, 1) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
567  answer.at(3, 2) = -( 2. * C1 * ( 6. * F.at(2, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 2) - 2. * I1 * invF.at(2, 2) * invF.at(3, 3) - 3. * I1 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
568  answer.at(3, 3) = ( 2. * C1 * ( 5. * I1 * invF.at(3, 3) * invF.at(3, 3) - 12. * F.at(3, 3) * invF.at(3, 3) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
569  answer.at(3, 4) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 3) - 3. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
570  answer.at(3, 5) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
571 
572  answer.at(4, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(2, 1) + 6. * F.at(1, 2) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
573  answer.at(4, 2) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(2, 1) - 5. * I1 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
574  answer.at(4, 3) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(2, 1) - 2. * I1 * invF.at(2, 1) * invF.at(3, 3) - 3. * I1 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 2. / 3.) );
575  answer.at(4, 4) = ( 2. * C1 * ( 5. * I1 * invF.at(2, 1) * invF.at(2, 1) - 12. * F.at(1, 2) * invF.at(2, 1) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
576  answer.at(4, 5) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
577 
578 
579  answer.at(5, 1) = -( 2. * C1 * ( 6. * F.at(1, 1) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(1, 1) - 5. * I1 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
580  answer.at(5, 2) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(2, 2) + 6. * F.at(2, 2) * invF.at(1, 2) - 5. * I1 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
581  answer.at(5, 3) = -( 2. * C1 * ( 6. * F.at(2, 1) * invF.at(3, 3) + 6. * F.at(3, 3) * invF.at(1, 2) - 2. * I1 * invF.at(1, 2) * invF.at(3, 3) - 3. * I1 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 2. / 3.) );
582  answer.at(5, 4) = -( 2. * C1 * ( 6. * F.at(1, 2) * invF.at(1, 2) + 6. * F.at(2, 1) * invF.at(2, 1) - 3. * I1 * invF.at(1, 1) * invF.at(2, 2) - 2. * I1 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 2. / 3.) );
583  answer.at(5, 5) = ( 2. * C1 * ( 5. * I1 * invF.at(1, 2) * invF.at(1, 2) - 12. * F.at(2, 1) * invF.at(1, 2) + 9 ) ) / ( 9. * pow(J, 2. / 3.) );
584 
585 
586 
588 
589 
590  answer.at(1, 1) = answer.at(1, 1) + ( 2. * C2 * ( 9. * F.at(1, 1) * F.at(1, 1) - 24. * I1 * F.at(1, 1) * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 1) + 24. * FC.at(1, 1) * invF.at(1, 1) + 9. * B.at(1, 1) - 9. * C.at(1, 1) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
591  answer.at(1, 2) = answer.at(1, 2) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 1) + 18. * F.at(1, 1) * F.at(2, 2) - 9. * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(1, 1) - 8. * I2 * invF.at(1, 1) * invF.at(2, 2) + 6. * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
592  answer.at(1, 3) = answer.at(1, 3) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 1) - 8. * I2 * invF.at(1, 1) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
593  answer.at(1, 4) = answer.at(1, 4) - ( 2. * C2 * ( 9. * C.at(2, 1) - 12. * FC.at(1, 1) * invF.at(1, 2) - 12. * FC.at(1, 2) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 2) + 12. * F.at(1, 1) * I1 * invF.at(2, 1) + 12. * F.at(1, 2) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
594  answer.at(1, 5) = answer.at(1, 5) + ( 2. * C2 * ( 9. * B.at(1, 2) + 12. * FC.at(1, 1) * invF.at(2, 1) + 12. * FC.at(2, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
595 
596  answer.at(2, 1) = answer.at(2, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(2, 2) - 9. * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(2, 2) + 6 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
597  answer.at(2, 2) = answer.at(2, 2) + ( 2. * C2 * ( 9. * F.at(2, 2) * F.at(2, 2) - 24. * I1 * F.at(2, 2) * invF.at(2, 2) - 2. * I2 * invF.at(2, 2) * invF.at(2, 2) + 24. * FC.at(2, 2) * invF.at(2, 2) + 9. * B.at(2, 2) - 9. * C.at(2, 2) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
598  answer.at(2, 3) = answer.at(2, 3) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 2) + 18 * F.at(2, 2) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 3) + 6 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
599  answer.at(2, 4) = answer.at(2, 4) + ( 2. * C2 * ( 9. * B.at(2, 1) + 12. * FC.at(1, 2) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(1, 2) + 9. * F.at(1, 2) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
600  answer.at(2, 5) = answer.at(2, 5) - ( 2. * C2 * ( 9. * C.at(1, 2) - 12. * FC.at(2, 1) * invF.at(2, 2) - 12. * FC.at(2, 2) * invF.at(2, 1) - 9. * F.at(2, 1) * F.at(2, 2) + 12. * F.at(2, 1) * I1 * invF.at(2, 2) + 12. * F.at(2, 2) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
601 
602  answer.at(3, 1) = answer.at(3, 1) + ( 2. * C2 * ( 12. * FC.at(1, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 1) + 18 * F.at(1, 1) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 1) - 12. * F.at(1, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 1) - 8 * I2 * invF.at(1, 1) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
603  answer.at(3, 2) = answer.at(3, 2) + ( 2. * C2 * ( 12. * FC.at(2, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 2) + 18 * F.at(2, 2) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 2) - 12. * F.at(2, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 2) - 8 * I2 * invF.at(2, 2) * invF.at(3, 3) + 6 * I2 * invF.at(2, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
604  answer.at(3, 3) = answer.at(3, 3) + ( 2. * C2 * ( 9. * F.at(3, 3) * F.at(3, 3) - 24. * I1 * F.at(3, 3) * invF.at(3, 3) - 2. * I2 * invF.at(3, 3) * invF.at(3, 3) + 24. * FC.at(3, 3) * invF.at(3, 3) + 9. * B.at(3, 3) - 9. * C.at(3, 3) + 9. * I1 ) ) / ( 9. * pow(J, 4. / 3.) );
605  answer.at(3, 4) = answer.at(3, 4) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 2) + 18 * F.at(1, 2) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 3) + 6 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
606  answer.at(3, 5) = answer.at(3, 5) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 1) + 18 * F.at(2, 1) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
607 
608 
609  answer.at(4, 1) = answer.at(4, 1) - ( 2. * C2 * ( 9. * C.at(1, 2) - 12. * FC.at(1, 1) * invF.at(2, 1) - 12. * FC.at(1, 2) * invF.at(1, 1) - 9. * F.at(1, 1) * F.at(1, 2) + 12. * F.at(1, 1) * I1 * invF.at(2, 1) + 12. * F.at(1, 2) * I1 * invF.at(1, 1) + 2. * I2 * invF.at(1, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
610  answer.at(4, 2) = answer.at(4, 2) + ( 2. * C2 * ( 9. * B.at(1, 2) + 12. * FC.at(1, 2) * invF.at(2, 2) + 12. * FC.at(2, 2) * invF.at(2, 1) + 9. * F.at(1, 2) * F.at(2, 2) - 12. * F.at(1, 2) * I1 * invF.at(2, 2) - 12. * F.at(2, 2) * I1 * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
611  answer.at(4, 3) = answer.at(4, 3) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(2, 1) + 18 * F.at(1, 2) * F.at(3, 3) - 9. * F.at(1, 3) * F.at(3, 2) - 12. * F.at(1, 2) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(2, 1) - 8 * I2 * invF.at(2, 1) * invF.at(3, 3) + 6 * I2 * invF.at(3, 1) * invF.at(2, 3) ) ) / ( 9. * pow(J, 4. / 3.) );
612  answer.at(4, 4) = answer.at(4, 4) + ( 2. * C2 * ( 9. * F.at(1, 2) * F.at(1, 2) - 24. * I1 * F.at(1, 2) * invF.at(2, 1) - 2. * I2 * invF.at(2, 1) * invF.at(2, 1) + 12. * FC.at(1, 2) * invF.at(2, 1) + 9. * B.at(1, 1) - 9. * C.at(2, 2) + 9. * I1 + 12. * FC.at(1, 2) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
613  answer.at(4, 5) = answer.at(4, 5) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(2, 1) + 12. * FC.at(2, 1) * invF.at(2, 1) - 9. * F.at(1, 1) * F.at(2, 2) + 18 * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
614 
615 
616 
617  answer.at(5, 1) = answer.at(5, 1) + ( 2. * C2 * ( 9. * B.at(2, 1) + 12. * FC.at(1, 1) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 1) + 9. * F.at(1, 1) * F.at(2, 1) - 12. * F.at(1, 1) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(1, 1) - 2. * I2 * invF.at(1, 1) * invF.at(1, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
618  answer.at(5, 2) = answer.at(5, 2) - ( 2. * C2 * ( 9. * C.at(2, 1) - 12. * FC.at(2, 1) * invF.at(2, 2) - 12. * FC.at(2, 2) * invF.at(1, 2) - 9. * F.at(2, 1) * F.at(2, 2) + 12. * F.at(2, 1) * I1 * invF.at(2, 2) + 12. * F.at(2, 2) * I1 * invF.at(1, 2) + 2. * I2 * invF.at(1, 2) * invF.at(2, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
619  answer.at(5, 3) = answer.at(5, 3) + ( 2. * C2 * ( 12. * FC.at(2, 1) * invF.at(3, 3) + 12. * FC.at(3, 3) * invF.at(1, 2) + 18 * F.at(2, 1) * F.at(3, 3) - 9. * F.at(2, 3) * F.at(3, 1) - 12. * F.at(2, 1) * I1 * invF.at(3, 3) - 12. * F.at(3, 3) * I1 * invF.at(1, 2) - 8 * I2 * invF.at(1, 2) * invF.at(3, 3) + 6 * I2 * invF.at(1, 3) * invF.at(3, 2) ) ) / ( 9. * pow(J, 4. / 3.) );
620  answer.at(5, 4) = answer.at(5, 4) + ( 2. * C2 * ( 12. * FC.at(1, 2) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 2) - 9. * F.at(1, 1) * F.at(2, 2) + 18 * F.at(1, 2) * F.at(2, 1) - 12. * F.at(1, 2) * I1 * invF.at(1, 2) - 12. * F.at(2, 1) * I1 * invF.at(2, 1) + 6 * I2 * invF.at(1, 1) * invF.at(2, 2) - 8 * I2 * invF.at(1, 2) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
621  answer.at(5, 5) = answer.at(5, 5) + ( 2. * C2 * ( 9. * F.at(2, 1) * F.at(2, 1) - 24. * I1 * F.at(2, 1) * invF.at(1, 2) - 2. * I2 * invF.at(1, 2) * invF.at(1, 2) + 12. * FC.at(2, 1) * invF.at(1, 2) + 9. * B.at(2, 2) - 9. * C.at(1, 1) + 9. * I1 + 12. * FC.at(2, 1) * invF.at(2, 1) ) ) / ( 9. * pow(J, 4. / 3.) );
622 
624  answer.at(1, 1) = answer.at(1, 1) - K *invF.at(1, 1) * invF.at(1, 1) * ( lnJ - 1. );
625  answer.at(1, 2) = answer.at(1, 2) + K * ( invF.at(1, 1) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 1) * lnJ );
626  answer.at(1, 3) = answer.at(1, 3) + K * ( invF.at(1, 1) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 1) * lnJ );
627  answer.at(1, 4) = answer.at(1, 4) - K *invF.at(1, 1) * invF.at(2, 1) * ( lnJ - 1. );
628  answer.at(1, 5) = answer.at(1, 5) - K *invF.at(1, 1) * invF.at(1, 2) * ( lnJ - 1. );
629 
630  answer.at(2, 1) = answer.at(2, 1) + K * ( invF.at(1, 1) * invF.at(2, 2) - invF.at(1, 2) * invF.at(2, 1) * lnJ );
631  answer.at(2, 2) = answer.at(2, 2) - K *invF.at(2, 2) * invF.at(2, 2) * ( lnJ - 1. );
632  answer.at(2, 3) = answer.at(2, 3) + K * ( invF.at(2, 2) * invF.at(3, 3) - invF.at(2, 3) * invF.at(3, 2) * lnJ );
633  answer.at(2, 4) = answer.at(2, 4) - K *invF.at(2, 1) * invF.at(2, 2) * ( lnJ - 1. );
634  answer.at(2, 5) = answer.at(2, 5) - K *invF.at(1, 2) * invF.at(2, 2) * ( lnJ - 1. );
635 
636  answer.at(3, 1) = answer.at(3, 1) + K * ( invF.at(1, 1) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 1) * lnJ );
637  answer.at(3, 2) = answer.at(3, 2) + K * ( invF.at(2, 2) * invF.at(3, 3) - invF.at(2, 3) * invF.at(3, 2) * lnJ );
638  answer.at(3, 3) = answer.at(3, 3) - K *invF.at(3, 3) * invF.at(3, 3) * ( lnJ - 1. );
639  answer.at(3, 4) = answer.at(3, 4) + K * ( invF.at(2, 1) * invF.at(3, 3) - invF.at(3, 1) * invF.at(2, 3) * lnJ );
640  answer.at(3, 5) = answer.at(3, 5) + K * ( invF.at(1, 2) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 2) * lnJ );
641 
642  answer.at(4, 1) = answer.at(4, 1) - K *invF.at(1, 1) * invF.at(2, 1) * ( lnJ - 1. );
643  answer.at(4, 2) = answer.at(4, 2) - K *invF.at(2, 1) * invF.at(2, 2) * ( lnJ - 1. );
644  answer.at(4, 3) = answer.at(4, 3) + K * ( invF.at(2, 1) * invF.at(3, 3) - invF.at(3, 1) * invF.at(2, 3) * lnJ );
645  answer.at(4, 4) = answer.at(4, 4) - K *invF.at(2, 1) * invF.at(2, 1) * ( lnJ - 1. );
646  answer.at(4, 5) = answer.at(4, 5) + K * ( invF.at(1, 2) * invF.at(2, 1) - invF.at(1, 1) * invF.at(2, 2) * lnJ );
647 
648  answer.at(5, 1) = answer.at(5, 1) - K *invF.at(1, 1) * invF.at(1, 2) * ( lnJ - 1. );
649  answer.at(5, 2) = answer.at(5, 2) - K *invF.at(1, 2) * invF.at(2, 2) * ( lnJ - 1. );
650  answer.at(5, 3) = answer.at(5, 3) + K * ( invF.at(1, 2) * invF.at(3, 3) - invF.at(1, 3) * invF.at(3, 2) * lnJ );
651  answer.at(5, 4) = answer.at(5, 4) + K * ( invF.at(1, 2) * invF.at(2, 1) - invF.at(1, 1) * invF.at(2, 2) * lnJ );
652  answer.at(5, 5) = answer.at(5, 5) - K *invF.at(1, 2) * invF.at(1, 2) * ( lnJ - 1. );
653 }
654 
655 
656 
657 MaterialStatus *
658 MooneyRivlinMaterial :: CreateStatus(GaussPoint *gp) const
659 {
660  return new StructuralMaterialStatus(1, this->giveDomain(), gp);
661 }
662 
663 
664 IRResultType
665 MooneyRivlinMaterial :: initializeFrom(InputRecord *ir)
666 {
667  IRResultType result; // Required by IR_GIVE_FIELD macro
668 
669 
670  result = StructuralMaterial :: initializeFrom(ir);
671  if ( result != IRRT_OK ) return result;
672 
673  IR_GIVE_FIELD(ir, K, _IFT_MooneyRivlinMaterial_k);
674  IR_GIVE_FIELD(ir, C1, _IFT_MooneyRivlinMaterial_c1);
675  IR_GIVE_FIELD(ir, C2, _IFT_MooneyRivlinMaterial_c2);
676 
677  return IRRT_OK;
678 }
679 
680 } // end namespace oofem
oofem::FloatMatrix::giveDeterminant
double giveDeterminant() const
Returns the trace (sum of diagonal components) of the receiver.
Definition: floatmatrix.C:1408
oofem::MooneyRivlinMaterial::giveFirstPKStressVector_3d
virtual void giveFirstPKStressVector_3d(FloatArray &answer, GaussPoint *gp, const FloatArray &vF, TimeStep *tStep)
Default implementation relies on giveRealStressVector for second Piola-Kirchoff stress.
Definition: mooneyrivlin.C:49
oofem::Material::giveStatus
virtual MaterialStatus * giveStatus(GaussPoint *gp) const
Returns material status of receiver in given integration point.
Definition: material.C:244
oofem::MooneyRivlinMaterial::MooneyRivlinMaterial
MooneyRivlinMaterial(int n, Domain *d)
Definition: mooneyrivlin.C:44
oofem::Domain
Class and object Domain.
Definition: domain.h:115
oofem::StructuralMaterialStatus::letTempFVectorBe
void letTempFVectorBe(const FloatArray &v)
Assigns tempFVector to given vector v.
Definition: structuralms.h:143
_IFT_MooneyRivlinMaterial_k
#define _IFT_MooneyRivlinMaterial_k
Definition: mooneyrivlin.h:46
oofem::FloatArray::at
double & at(int i)
Coefficient access function.
Definition: floatarray.h:131
oofem::StructuralMaterialStatus
This class implements a structural material status information.
Definition: structuralms.h:65
_IFT_MooneyRivlinMaterial_c2
#define _IFT_MooneyRivlinMaterial_c2
Definition: mooneyrivlin.h:45
oofem::FloatArray::beVectorForm
void beVectorForm(const FloatMatrix &aMatrix)
Reciever will be a vector with 9 components formed from a 3x3 matrix.
Definition: floatarray.C:992
_IFT_MooneyRivlinMaterial_c1
#define _IFT_MooneyRivlinMaterial_c1
Definition: mooneyrivlin.h:44
oofem::MooneyRivlinMaterial::givePlaneStrainStiffMtrx_dPdF
virtual void givePlaneStrainStiffMtrx_dPdF(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Definition: mooneyrivlin.C:527
oofem::MatResponseMode
MatResponseMode
Describes the character of characteristic material matrix.
Definition: matresponsemode.h:64
oofem::MooneyRivlinMaterial::CreateStatus
virtual MaterialStatus * CreateStatus(GaussPoint *gp) const
Creates new copy of associated status and inserts it into given integration point.
Definition: mooneyrivlin.C:658
oofem::MooneyRivlinMaterial::giveFirstPKStressVector_PlaneStrain
void giveFirstPKStressVector_PlaneStrain(FloatArray &answer, GaussPoint *gp, const FloatArray &reducedvF, TimeStep *tStep)
Default implementation relies on giveFirstPKStressVector_3d.
Definition: mooneyrivlin.C:96
oofem::FloatMatrix::beMatrixForm
void beMatrixForm(const FloatArray &aArray)
Definition: floatmatrix.C:1657
oofem::StructuralMaterial::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: structuralmaterial.C:2554
oofem::StructuralMaterial::giveReducedVectorForm
static void giveReducedVectorForm(FloatArray &answer, const FloatArray &vec, MaterialMode matMode)
Converts the full symmetric Voigt vector (2nd order tensor) to reduced form.
Definition: structuralmaterial.C:2488
oofem::FloatMatrix::at
double at(int i, int j) const
Coefficient access function.
Definition: floatmatrix.h:176
oofem::MooneyRivlinMaterial::give3dMaterialStiffnessMatrix_dPdF
virtual void give3dMaterialStiffnessMatrix_dPdF(FloatMatrix &answer, MatResponseMode mode, GaussPoint *gp, TimeStep *tStep)
Definition: mooneyrivlin.C:148
oofem::MaterialStatus
Abstract base class representing a material status information.
Definition: matstatus.h:84
oofem::MooneyRivlinMaterial::initializeFrom
virtual IRResultType initializeFrom(InputRecord *ir)
Initializes receiver according to object description stored in input record.
Definition: mooneyrivlin.C:665
oofem::FloatArray
Class representing vector of real numbers.
Definition: floatarray.h:82
oofem::FloatMatrix
Implementation of matrix containing floating point numbers.
Definition: floatmatrix.h:94
oofem::IRResultType
IRResultType
Type defining the return values of InputRecord reading operations.
Definition: irresulttype.h:47
oofem::FloatMatrix::resize
void resize(int rows, int cols)
Checks size of receiver towards requested bounds.
Definition: floatmatrix.C:1358
oofem::InputRecord
Class representing the general Input Record.
Definition: inputrecord.h:101
oofem::FloatMatrix::add
void add(const FloatMatrix &a)
Adds matrix to the receiver.
Definition: floatmatrix.C:1023
oofem::FloatMatrix::beTProductOf
void beTProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:367
oofem::StructuralMaterialStatus::giveTempFVector
const FloatArray & giveTempFVector() const
Returns the const pointer to receiver&#39;s temporary deformation gradient vector.
Definition: structuralms.h:123
oofem::FloatMatrix::beTranspositionOf
void beTranspositionOf(const FloatMatrix &src)
Assigns to the receiver the transposition of parameter.
Definition: floatmatrix.C:323
oofem::StructuralMaterial
Abstract base class for all "structural" constitutive models.
Definition: structuralmaterial.h:113
oofem::FloatMatrix::beProductTOf
void beProductTOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:397
oofem::FEMComponent::giveDomain
Domain * giveDomain() const
Definition: femcmpnn.h:100
oofem::REGISTER_Material
REGISTER_Material(DummyMaterial)
oofem::FloatMatrix::beProductOf
void beProductOf(const FloatMatrix &a, const FloatMatrix &b)
Assigns to the receiver product of .
Definition: floatmatrix.C:337
oofem
the oofem namespace is to define a context or scope in which all oofem names are defined.
IR_GIVE_FIELD
#define IR_GIVE_FIELD(__ir, __value, __id)
Macro facilitating the use of input record reading methods.
Definition: inputrecord.h:69
oofem::FloatMatrix::beInverseOf
void beInverseOf(const FloatMatrix &src)
Modifies receiver to become inverse of given parameter.
Definition: floatmatrix.C:835
oofem::GaussPoint
Class representing integration point in finite element program.
Definition: gausspoint.h:93
oofem::TimeStep
Class representing solution step.
Definition: timestep.h:80
oofem::StructuralMaterialStatus::letTempPVectorBe
void letTempPVectorBe(const FloatArray &v)
Assigns tempPVector to given vector v.
Definition: structuralms.h:139
This page is part of the OOFEM documentation. Copyright (c) 2011 Borek Patzak
Project e-mail: info@oofem.org
Generated at Tue Jan 2 2018 20:07:30 for OOFEM by doxygen 1.8.11 written by Dimitri van Heesch, © 1997-2011