[cig-commits] r16429 - in short/3D/PyLith/trunk/libsrc: feassemble materials
willic3 at geodynamics.org
willic3 at geodynamics.org
Tue Mar 16 18:03:37 PDT 2010
Author: willic3
Date: 2010-03-16 18:03:36 -0700 (Tue, 16 Mar 2010)
New Revision: 16429
Added:
short/3D/PyLith/trunk/libsrc/feassemble/jacobian2d_nonsymm_lgdeform.wxm
short/3D/PyLith/trunk/libsrc/feassemble/jacobian3d_nonsymm_lgdeform.wxm
Modified:
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc
short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.hh
short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.icc
short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
Log:
Changed everything to handle full elasticity matrix, including the
possibility of a nonsymmetric matrix.
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticity.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -817,7 +817,7 @@
{ // _elasticityJacobian2D
const int spaceDim = 2;
const int cellDim = 2;
- const int numConsts = 6;
+ const int numConsts = 9;
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -838,9 +838,12 @@
const double C1111 = elasticConsts[iQuad*numConsts+0];
const double C1122 = elasticConsts[iQuad*numConsts+1];
const double C1112 = elasticConsts[iQuad*numConsts+2]/2.0;
- const double C2222 = elasticConsts[iQuad*numConsts+3];
- const double C2212 = elasticConsts[iQuad*numConsts+4]/2.0;
- const double C1212 = elasticConsts[iQuad*numConsts+5]/2.0;
+ const double C2211 = elasticConsts[iQuad*numConsts+3];
+ const double C2222 = elasticConsts[iQuad*numConsts+4];
+ const double C2212 = elasticConsts[iQuad*numConsts+5]/2.0;
+ const double C1211 = elasticConsts[iQuad*numConsts+6]/2.0;
+ const double C1222 = elasticConsts[iQuad*numConsts+7]/2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+8]/2.0;
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
@@ -852,16 +855,16 @@
const double Njp = basisDeriv[iQ+jBasis*spaceDim ];
const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp +
+ C1111 * Nip * Njp + C1211 * Niq * Njp +
C1112 * Nip * Njq + C1212 * Niq * Njq;
const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq +
+ C1122 * Nip * Njq + C1222 * Niq * Njq +
C1112 * Nip * Njp + C1212 * Niq * Njp;
const double ki1j0 =
- C1122 * Niq * Njp + C2212 * Niq * Njq +
- C1112 * Nip * Njp + C1212 * Nip * Njq;
+ C2211 * Niq * Njp + C1211 * Nip * Njp +
+ C2212 * Niq * Njq + C1212 * Nip * Njq;
const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq +
+ C2222 * Niq * Njq + C1222 * Nip * Njq +
C2212 * Niq * Njp + C1212 * Nip * Njp;
const int jBlock = (jBasis*spaceDim );
const int jBlock1 = (jBasis*spaceDim+1);
@@ -891,7 +894,6 @@
assert(2 == cellDim);
assert(quadWts.size() == numQuadPts);
- const int numConsts = 6;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const double wt = quadWts[iQuad] * jacobianDet[iQuad];
@@ -926,7 +928,7 @@
{ // _elasticityJacobian3D
const int spaceDim = 3;
const int cellDim = 3;
- const int numConsts = 21;
+ const int numConsts = 36;
const int numQuadPts = _quadrature->numQuadPts();
const int numBasis = _quadrature->numBasis();
@@ -951,21 +953,36 @@
const double C1112 = elasticConsts[iQuad*numConsts+ 3] / 2.0;
const double C1123 = elasticConsts[iQuad*numConsts+ 4] / 2.0;
const double C1113 = elasticConsts[iQuad*numConsts+ 5] / 2.0;
- const double C2222 = elasticConsts[iQuad*numConsts+ 6];
- const double C2233 = elasticConsts[iQuad*numConsts+ 7];
- const double C2212 = elasticConsts[iQuad*numConsts+ 8] / 2.0;
- const double C2223 = elasticConsts[iQuad*numConsts+ 9] / 2.0;
- const double C2213 = elasticConsts[iQuad*numConsts+10] / 2.0;
- const double C3333 = elasticConsts[iQuad*numConsts+11];
- const double C3312 = elasticConsts[iQuad*numConsts+12] / 2.0;
- const double C3323 = elasticConsts[iQuad*numConsts+13] / 2.0;
- const double C3313 = elasticConsts[iQuad*numConsts+14] / 2.0;
- const double C1212 = elasticConsts[iQuad*numConsts+15] / 2.0;
- const double C1223 = elasticConsts[iQuad*numConsts+16] / 2.0;
- const double C1213 = elasticConsts[iQuad*numConsts+17] / 2.0;
- const double C2323 = elasticConsts[iQuad*numConsts+18] / 2.0;
- const double C2313 = elasticConsts[iQuad*numConsts+19] / 2.0;
- const double C1313 = elasticConsts[iQuad*numConsts+20] / 2.0;
+ const double C2211 = elasticConsts[iQuad*numConsts+ 6];
+ const double C2222 = elasticConsts[iQuad*numConsts+ 7];
+ const double C2233 = elasticConsts[iQuad*numConsts+ 8];
+ const double C2212 = elasticConsts[iQuad*numConsts+ 9] / 2.0;
+ const double C2223 = elasticConsts[iQuad*numConsts+10] / 2.0;
+ const double C2213 = elasticConsts[iQuad*numConsts+11] / 2.0;
+ const double C3311 = elasticConsts[iQuad*numConsts+12];
+ const double C3322 = elasticConsts[iQuad*numConsts+13];
+ const double C3333 = elasticConsts[iQuad*numConsts+14];
+ const double C3312 = elasticConsts[iQuad*numConsts+15] / 2.0;
+ const double C3323 = elasticConsts[iQuad*numConsts+16] / 2.0;
+ const double C3313 = elasticConsts[iQuad*numConsts+17] / 2.0;
+ const double C1211 = elasticConsts[iQuad*numConsts+18] / 2.0;
+ const double C1222 = elasticConsts[iQuad*numConsts+19] / 2.0;
+ const double C1233 = elasticConsts[iQuad*numConsts+20] / 2.0;
+ const double C1212 = elasticConsts[iQuad*numConsts+21] / 2.0;
+ const double C1223 = elasticConsts[iQuad*numConsts+22] / 2.0;
+ const double C1213 = elasticConsts[iQuad*numConsts+23] / 2.0;
+ const double C2311 = elasticConsts[iQuad*numConsts+24] / 2.0;
+ const double C2322 = elasticConsts[iQuad*numConsts+25] / 2.0;
+ const double C2333 = elasticConsts[iQuad*numConsts+26] / 2.0;
+ const double C2312 = elasticConsts[iQuad*numConsts+27] / 2.0;
+ const double C2323 = elasticConsts[iQuad*numConsts+28] / 2.0;
+ const double C2313 = elasticConsts[iQuad*numConsts+29] / 2.0;
+ const double C1311 = elasticConsts[iQuad*numConsts+30] / 2.0;
+ const double C1322 = elasticConsts[iQuad*numConsts+31] / 2.0;
+ const double C1333 = elasticConsts[iQuad*numConsts+32] / 2.0;
+ const double C1312 = elasticConsts[iQuad*numConsts+33] / 2.0;
+ const double C1323 = elasticConsts[iQuad*numConsts+34] / 2.0;
+ const double C1313 = elasticConsts[iQuad*numConsts+35] / 2.0;
for (int iBasis=0, iQ=iQuad*numBasis*spaceDim;
iBasis < numBasis;
++iBasis) {
@@ -977,40 +994,40 @@
const double Njq = basisDeriv[iQ+jBasis*spaceDim+1];
const double Njr = basisDeriv[iQ+jBasis*spaceDim+2];
const double ki0j0 =
- C1111 * Nip * Njp + C1112 * Niq * Njp + C1113 * Nir * Njp +
- C1112 * Nip * Njq + C1212 * Niq * Njq + C1213 * Nir * Njq +
+ C1111 * Nip * Njp + C1211 * Niq * Njp + C1311 * Nir * Njp +
+ C1112 * Nip * Njq + C1212 * Niq * Njq + C1312 * Nir * Njq +
C1113 * Nip * Njr + C1213 * Niq * Njr + C1313 * Nir * Njr;
const double ki0j1 =
- C1122 * Nip * Njq + C2212 * Niq * Njq + C2213 * Nir * Njq +
- C1112 * Nip * Njp + C1212 * Niq * Njp + C1213 * Nir * Njp +
- C1123 * Nip * Njr + C1223 * Niq * Njr + C2313 * Nir * Njr;
+ C1122 * Nip * Njq + C1222 * Niq * Njq + C1322 * Nir * Njq +
+ C1112 * Nip * Njp + C1212 * Niq * Njp + C1312 * Nir * Njp +
+ C1123 * Nip * Njr + C1223 * Niq * Njr + C1323 * Nir * Njr;
const double ki0j2 =
- C1133 * Nip * Njr + C3312 * Niq * Njr + C3313 * Nir * Njr +
- C1123 * Nip * Njq + C1223 * Niq * Njq + C2313 * Nir * Njq +
+ C1133 * Nip * Njr + C1233 * Niq * Njr + C1333 * Nir * Njr +
+ C1123 * Nip * Njq + C1223 * Niq * Njq + C1323 * Nir * Njq +
C1113 * Nip * Njp + C1213 * Niq * Njp + C1313 * Nir * Njp;
const double ki1j0 =
- C1122 * Niq * Njp + C1112 * Nip * Njp + C1123 * Nir * Njp +
- C2212 * Niq * Njq + C1212 * Nip * Njq + C1223 * Nir * Njq +
+ C2211 * Niq * Njp + C1211 * Nip * Njp + C2311 * Nir * Njp +
+ C2212 * Niq * Njq + C1212 * Nip * Njq + C2312 * Nir * Njq +
C2213 * Niq * Njr + C1213 * Nip * Njr + C2313 * Nir * Njr;
const double ki1j1 =
- C2222 * Niq * Njq + C2212 * Nip * Njq + C2223 * Nir * Njq +
- C2212 * Niq * Njp + C1212 * Nip * Njp + C1223 * Nir * Njp +
+ C2222 * Niq * Njq + C1222 * Nip * Njq + C2322 * Nir * Njq +
+ C2212 * Niq * Njp + C1212 * Nip * Njp + C2312 * Nir * Njp +
C2223 * Niq * Njr + C1223 * Nip * Njr + C2323 * Nir * Njr;
const double ki1j2 =
- C2233 * Niq * Njr + C3312 * Nip * Njr + C3323 * Nir * Njr +
+ C2233 * Niq * Njr + C1233 * Nip * Njr + C2333 * Nir * Njr +
C2223 * Niq * Njq + C1223 * Nip * Njq + C2323 * Nir * Njq +
C2213 * Niq * Njp + C1213 * Nip * Njp + C2313 * Nir * Njp;
const double ki2j0 =
- C1133 * Nir * Njp + C1123 * Niq * Njp + C1113 * Nip * Njp +
- C3312 * Nir * Njq + C1223 * Niq * Njq + C1213 * Nip * Njq +
+ C3311 * Nir * Njp + C2311 * Niq * Njp + C1311 * Nip * Njp +
+ C3312 * Nir * Njq + C2312 * Niq * Njq + C1312 * Nip * Njq +
C3313 * Nir * Njr + C2313 * Niq * Njr + C1313 * Nip * Njr;
const double ki2j1 =
- C2233 * Nir * Njq + C2223 * Niq * Njq + C2213 * Nip * Njq +
- C3312 * Nir * Njp + C1223 * Niq * Njp + C1213 * Nip * Njp +
- C3323 * Nir * Njr + C2323 * Niq * Njr + C2313 * Nip * Njr;
+ C3322 * Nir * Njq + C2322 * Niq * Njq + C1322 * Nip * Njq +
+ C3312 * Nir * Njp + C2312 * Niq * Njp + C1312 * Nip * Njp +
+ C3323 * Nir * Njr + C2323 * Niq * Njr + C1323 * Nip * Njr;
const double ki2j2 =
- C3333 * Nir * Njr + C3323 * Niq * Njr + C3313 * Nip * Njr +
- C3323 * Nir * Njq + C2323 * Niq * Njq + C2313 * Nip * Njq +
+ C3333 * Nir * Njr + C2333 * Niq * Njr + C1333 * Nip * Njr +
+ C3323 * Nir * Njq + C2323 * Niq * Njq + C1323 * Nip * Njq +
C3313 * Nir * Njp + C2313 * Niq * Njp + C1313 * Nip * Njp;
const int iBlock = iBasis*spaceDim * (numBasis*spaceDim);
const int iBlock1 = (iBasis*spaceDim+1) * (numBasis*spaceDim);
@@ -1049,7 +1066,6 @@
assert(3 == cellDim);
assert(quadWts.size() == numQuadPts);
- const int numConsts = 21;
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
Modified: short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/feassemble/IntegratorElasticityLgDeform.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -559,7 +559,7 @@
assert(2 == cellDim);
assert(quadWts.size() == numQuadPts);
- const int numConsts = 6;
+ const int numConsts = 9;
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
const int iQ = iQuad*numBasis*spaceDim;
@@ -572,9 +572,12 @@
const double C1111 = elasticConsts[iC+0];
const double C1122 = elasticConsts[iC+1];
const double C1112 = elasticConsts[iC+2] / 2.0; // 2*mu -> mu
- const double C2222 = elasticConsts[iC+3];
- const double C2212 = elasticConsts[iC+4] / 2.0;
- const double C1212 = elasticConsts[iC+5] / 2.0;
+ const double C2211 = elasticConsts[iC+3];
+ const double C2222 = elasticConsts[iC+4];
+ const double C2212 = elasticConsts[iC+5] / 2.0;
+ const double C1211 = elasticConsts[iC+6] / 2.0;
+ const double C1222 = elasticConsts[iC+7] / 2.0;
+ const double C1212 = elasticConsts[iC+8] / 2.0;
const int iS = iQuad*tensorSize;
const double s11 = stress[iS+0];
@@ -613,40 +616,40 @@
const double Ki0j0 =
l12*Niq*(l12*Njq*C2222 +
((l11+1)*Njq+l12*Njp)*C2212 +
- (l11+1)*Njp*C1122) +
- ((l11+1)*Niq+l12*Nip)*(l12*Njq*C2212 +
+ (l11+1)*Njp*C2211) +
+ ((l11+1)*Niq+l12*Nip)*(l12*Njq*C1222 +
((l11+1)*Njq+l12*Njp)*C1212 +
- (l11+1)*Njp*C1112) +
+ (l11+1)*Njp*C1211) +
(l11+1)*Nip*(l12*Njq*C1122 +
((l11+1)*Njq+l12*Njp)*C1112 +
(l11+1)*Njp*C1111);
const double Ki0j1 =
l12*Niq*((l22+1.0)*Njq*C2222 +
(l21*Njq+(l22+1.0)*Njp)*C2212 +
- l21*Njp*C1122) +
- ((l11+1.0)*Niq+l12*Nip)*((l22+1.0)*Njq*C2212 +
+ l21*Njp*C2211) +
+ ((l11+1.0)*Niq+l12*Nip)*((l22+1.0)*Njq*C1222 +
(l21*Njq+(l22+1.0)*Njp)*C1212 +
- l21*Njp*C1112) +
+ l21*Njp*C1211) +
(l11+1.0)*Nip*((l22+1.0)*Njq*C1122 +
(l21*Njq+(l22+1.0)*Njp)*C1112 +
l21*Njp*C1111);
const double Ki1j0 =
(l22+1.0)*Niq*(l12*Njq*C2222 +
((l11+1.0)*Njq+l12*Njp)*C2212 +
- (l11+1.0)*Njp*C1122) +
- (l21*Niq+(l22+1.0)*Nip)*(l12*Njq*C2212 +
+ (l11+1.0)*Njp*C2211) +
+ (l21*Niq+(l22+1.0)*Nip)*(l12*Njq*C1222 +
((l11+1.0)*Njq+l12*Njp)*C1212 +
- (l11+1.0)*Njp*C1112) +
+ (l11+1.0)*Njp*C1211) +
l21*Nip*(l12*Njq*C1122 +
((l11+1.0)*Njq+l12*Njp)*C1112 +
(l11+1.0)*Njp*C1111);
const double Ki1j1 =
(l22+1.0)*Niq*((l22+1.0)*Njq*C2222 +
(l21*Njq+(l22+1.0)*Njp)*C2212 +
- l21*Njp*C1122) +
- (l21*Niq+(l22+1.0)*Nip)*((l22+1.0)*Njq*C2212 +
+ l21*Njp*C2211) +
+ (l21*Niq+(l22+1.0)*Nip)*((l22+1.0)*Njq*C1222 +
(l21*Njq+(l22+1.0)*Njp)*C1212 +
- l21*Njp*C1112) +
+ l21*Njp*C1211) +
l21*Nip*((l22+1.0)*Njq*C1122 +
(l21*Njq+(l22+1.0)*Njp)*C1112 +
l21*Njp*C1111);
@@ -685,7 +688,7 @@
assert(3 == cellDim);
assert(quadWts.size() == numQuadPts);
assert(6 == tensorSize);
- const int numConsts = 21;
+ const int numConsts = 36;
// Compute Jacobian for consistent tangent matrix
for (int iQuad=0; iQuad < numQuadPts; ++iQuad) {
@@ -702,21 +705,36 @@
const double C1112 = elasticConsts[iC+ 3] / 2.0;
const double C1123 = elasticConsts[iC+ 4] / 2.0;
const double C1113 = elasticConsts[iC+ 5] / 2.0;
- const double C2222 = elasticConsts[iC+ 6];
- const double C2233 = elasticConsts[iC+ 7];
- const double C2212 = elasticConsts[iC+ 8] / 2.0;
- const double C2223 = elasticConsts[iC+ 9] / 2.0;
- const double C2213 = elasticConsts[iC+10] / 2.0;
- const double C3333 = elasticConsts[iC+11];
- const double C3312 = elasticConsts[iC+12] / 2.0;
- const double C3323 = elasticConsts[iC+13] / 2.0;
- const double C3313 = elasticConsts[iC+14] / 2.0;
- const double C1212 = elasticConsts[iC+15] / 2.0;
- const double C1223 = elasticConsts[iC+16] / 2.0;
- const double C1213 = elasticConsts[iC+17] / 2.0;
- const double C2323 = elasticConsts[iC+18] / 2.0;
- const double C2313 = elasticConsts[iC+19] / 2.0;
- const double C1313 = elasticConsts[iC+20] / 2.0;
+ const double C2211 = elasticConsts[iC+ 6];
+ const double C2222 = elasticConsts[iC+ 7];
+ const double C2233 = elasticConsts[iC+ 8];
+ const double C2212 = elasticConsts[iC+ 9] / 2.0;
+ const double C2223 = elasticConsts[iC+10] / 2.0;
+ const double C2213 = elasticConsts[iC+11] / 2.0;
+ const double C3311 = elasticConsts[iC+12];
+ const double C3322 = elasticConsts[iC+13];
+ const double C3333 = elasticConsts[iC+14];
+ const double C3312 = elasticConsts[iC+15] / 2.0;
+ const double C3323 = elasticConsts[iC+16] / 2.0;
+ const double C3313 = elasticConsts[iC+17] / 2.0;
+ const double C1211 = elasticConsts[iC+18] / 2.0;
+ const double C1222 = elasticConsts[iC+19] / 2.0;
+ const double C1233 = elasticConsts[iC+20] / 2.0;
+ const double C1212 = elasticConsts[iC+21] / 2.0;
+ const double C1223 = elasticConsts[iC+22] / 2.0;
+ const double C1213 = elasticConsts[iC+23] / 2.0;
+ const double C2311 = elasticConsts[iC+24] / 2.0;
+ const double C2322 = elasticConsts[iC+25] / 2.0;
+ const double C2333 = elasticConsts[iC+26] / 2.0;
+ const double C2312 = elasticConsts[iC+27] / 2.0;
+ const double C2323 = elasticConsts[iC+28] / 2.0;
+ const double C2313 = elasticConsts[iC+29] / 2.0;
+ const double C1311 = elasticConsts[iC+30] / 2.0;
+ const double C1322 = elasticConsts[iC+31] / 2.0;
+ const double C1333 = elasticConsts[iC+32] / 2.0;
+ const double C1312 = elasticConsts[iC+33] / 2.0;
+ const double C1323 = elasticConsts[iC+34] / 2.0;
+ const double C1313 = elasticConsts[iC+35] / 2.0;
const int iS = iQuad*tensorSize;
const double s11 = stress[iS+0];
@@ -767,32 +785,32 @@
(l12*Njr+l13*Njq)*C3323 +
((l11+1)*Njr+l13*Njp)*C3313 +
((l11+1)*Njq+l12*Njp)*C3312 +
- l12*Njq*C2233 +
- (l11+1)*Njp*C1133) +
- (l12*Nir+l13*Niq)*(l13*Njr*C3323 +
+ l12*Njq*C3322 +
+ (l11+1)*Njp*C3311) +
+ (l12*Nir+l13*Niq)*(l13*Njr*C2333 +
(l12*Njr+l13*Njq)*C2323 +
((l11+1)*Njr+l13*Njp)*C2313 +
- l12*Njq*C2223 +
- ((l11+1)*Njq+l12*Njp)*C1223 +
- (l11+1)*Njp*C1123) +
- ((l11+1)*Nir+l13*Nip)*(l13*Njr*C3313 +
- (l12*Njr+l13*Njq)*C2313 +
- l12*Njq*C2213 +
+ l12*Njq*C2322 +
+ ((l11+1)*Njq+l12*Njp)*C2312 +
+ (l11+1)*Njp*C2311) +
+ ((l11+1)*Nir+l13*Nip)*(l13*Njr*C1333 +
+ (l12*Njr+l13*Njq)*C1323 +
+ l12*Njq*C1322 +
((l11+1)*Njr+l13*Njp)*C1313 +
- ((l11+1)*Njq+l12*Njp)*C1213 +
- (l11+1)*Njp*C1113) +
- ((l11+1)*Niq+l12*Nip)*(l13*Njr*C3312 +
- l12*Njq*C2212 +
+ ((l11+1)*Njq+l12*Njp)*C1312 +
+ (l11+1)*Njp*C1311) +
+ ((l11+1)*Niq+l12*Nip)*(l13*Njr*C1233 +
+ l12*Njq*C1222 +
(l12*Njr+l13*Njq)*C1223 +
((l11+1)*Njr+l13*Njp)*C1213 +
((l11+1)*Njq+l12*Njp)*C1212 +
- (l11+1)*Njp*C1112) +
+ (l11+1)*Njp*C1211) +
l12*Niq*(l13*Njr*C2233 +
(l12*Njr+l13*Njq)*C2223 +
l12*Njq*C2222 +
((l11+1)*Njr+l13*Njp)*C2213 +
((l11+1)*Njq+l12*Njp)*C2212 +
- (l11+1)*Njp*C1122) +
+ (l11+1)*Njp*C2211) +
(l11+1)*Nip*(l13*Njr*C1133 +
(l12*Njr+l13*Njq)*C1123 +
l12*Njq*C1122 +
@@ -805,32 +823,32 @@
((l22+1)*Njr+l23*Njq)*C3323 +
(l21*Njr+l23*Njp)*C3313 +
(l21*Njq+(l22+1)*Njp)*C3312 +
- (l22+1)*Njq*C2233 +
- l21*Njp*C1133) +
- (l12*Nir+l13*Niq)*(l23*Njr*C3323 +
+ (l22+1)*Njq*C3322 +
+ l21*Njp*C3311) +
+ (l12*Nir+l13*Niq)*(l23*Njr*C2333 +
((l22+1)*Njr+l23*Njq)*C2323 +
(l21*Njr+l23*Njp)*C2313 +
- (l22+1)*Njq*C2223 +
- (l21*Njq+(l22+1)*Njp)*C1223 +
- l21*Njp*C1123) +
- ((l11+1)*Nir+l13*Nip)*(l23*Njr*C3313 +
- ((l22+1)*Njr+l23*Njq)*C2313 +
- (l22+1)*Njq*C2213 +
+ (l22+1)*Njq*C2322 +
+ (l21*Njq+(l22+1)*Njp)*C2312 +
+ l21*Njp*C2311) +
+ ((l11+1)*Nir+l13*Nip)*(l23*Njr*C1333 +
+ ((l22+1)*Njr+l23*Njq)*C1323 +
+ (l22+1)*Njq*C1322 +
(l21*Njr+l23*Njp)*C1313 +
- (l21*Njq+(l22+1)*Njp)*C1213 +
- l21*Njp*C1113) +
- ((l11+1)*Niq+l12*Nip)*(l23*Njr*C3312 +
- (l22+1)*Njq*C2212 +
+ (l21*Njq+(l22+1)*Njp)*C1312 +
+ l21*Njp*C1311) +
+ ((l11+1)*Niq+l12*Nip)*(l23*Njr*C1233 +
+ (l22+1)*Njq*C1222 +
((l22+1)*Njr+l23*Njq)*C1223 +
(l21*Njr+l23*Njp)*C1213 +
(l21*Njq+(l22+1)*Njp)*C1212 +
- l21*Njp*C1112) +
+ l21*Njp*C1211) +
l12*Niq*(l23*Njr*C2233 +
((l22+1)*Njr+l23*Njq)*C2223 +
(l22+1)*Njq*C2222 +
(l21*Njr+l23*Njp)*C2213 +
(l21*Njq+(l22+1)*Njp)*C2212 +
- l21*Njp*C1122) +
+ l21*Njp*C2211) +
(l11+1)*Nip*(l23*Njr*C1133 +
((l22+1)*Njr+l23*Njq)*C1123 +
(l22+1)*Njq*C1122 +
@@ -843,32 +861,32 @@
(l32*Njr+(l33+1)*Njq)*C3323 +
(l31*Njr+(l33+1)*Njp)*C3313 +
(l31*Njq+l32*Njp)*C3312 +
- l32*Njq*C2233 +
- l31*Njp*C1133) +
- (l12*Nir+l13*Niq)*((l33+1)*Njr*C3323 +
+ l32*Njq*C3322 +
+ l31*Njp*C3311) +
+ (l12*Nir+l13*Niq)*((l33+1)*Njr*C2333 +
(l32*Njr+(l33+1)*Njq)*C2323 +
(l31*Njr+(l33+1)*Njp)*C2313 +
- l32*Njq*C2223 +
- (l31*Njq+l32*Njp)*C1223 +
- l31*Njp*C1123) +
- ((l11+1)*Nir+l13*Nip)*((l33+1)*Njr*C3313 +
- (l32*Njr+(l33+1)*Njq)*C2313 +
- l32*Njq*C2213 +
+ l32*Njq*C2322 +
+ (l31*Njq+l32*Njp)*C2312 +
+ l31*Njp*C2311) +
+ ((l11+1)*Nir+l13*Nip)*((l33+1)*Njr*C1333 +
+ (l32*Njr+(l33+1)*Njq)*C1323 +
+ l32*Njq*C1322 +
(l31*Njr+(l33+1)*Njp)*C1313 +
- (l31*Njq+l32*Njp)*C1213 +
- l31*Njp*C1113) +
- ((l11+1)*Niq+l12*Nip)*((l33+1)*Njr*C3312 +
- l32*Njq*C2212 +
+ (l31*Njq+l32*Njp)*C1312 +
+ l31*Njp*C1311) +
+ ((l11+1)*Niq+l12*Nip)*((l33+1)*Njr*C1233 +
+ l32*Njq*C1222 +
(l32*Njr+(l33+1)*Njq)*C1223 +
(l31*Njr+(l33+1)*Njp)*C1213 +
(l31*Njq+l32*Njp)*C1212 +
- l31*Njp*C1112) +
+ l31*Njp*C1211) +
l12*Niq*((l33+1)*Njr*C2233 +
(l32*Njr+(l33+1)*Njq)*C2223 +
l32*Njq*C2222 +
(l31*Njr+(l33+1)*Njp)*C2213 +
(l31*Njq+l32*Njp)*C2212 +
- l31*Njp*C1122) +
+ l31*Njp*C2211) +
(l11+1)*Nip*((l33+1)*Njr*C1133 +
(l32*Njr+(l33+1)*Njq)*C1123 +
l32*Njq*C1122 +
@@ -881,32 +899,32 @@
(l12*Njr+l13*Njq)*C3323 +
((l11+1)*Njr+l13*Njp)*C3313 +
((l11+1)*Njq+l12*Njp)*C3312 +
- l12*Njq*C2233+(l11+1)*Njp*C1133) +
- ((l22+1)*Nir+l23*Niq)*(l13*Njr*C3323 +
+ l12*Njq*C3322 +
+ (l11+1)*Njp*C3311) +
+ ((l22+1)*Nir+l23*Niq)*(l13*Njr*C2333 +
(l12*Njr+l13*Njq)*C2323 +
((l11+1)*Njr+l13*Njp)*C2313 +
- l12*Njq*C2223 +
- ((l11+1)*Njq+l12*Njp)*C1223 +
- (l11+1)*Njp*C1123) +
- (l21*Nir+l23*Nip)*(l13*Njr*C3313 +
- (l12*Njr+l13*Njq)*C2313 +
- l12*Njq*C2213 +
- ((l11+1)*Njr +
- l13*Njp)*C1313 +
- ((l11+1)*Njq+l12*Njp)*C1213 +
- (l11+1)*Njp*C1113) +
- (l21*Niq+(l22+1)*Nip)*(l13*Njr*C3312 +
- l12*Njq*C2212 +
+ l12*Njq*C2322 +
+ ((l11+1)*Njq+l12*Njp)*C2312 +
+ (l11+1)*Njp*C2311) +
+ (l21*Nir+l23*Nip)*(l13*Njr*C1333 +
+ (l12*Njr+l13*Njq)*C1323 +
+ l12*Njq*C1322 +
+ ((l11+1)*Njr+l13*Njp)*C1313 +
+ ((l11+1)*Njq+l12*Njp)*C1312 +
+ (l11+1)*Njp*C1311) +
+ (l21*Niq+(l22+1)*Nip)*(l13*Njr*C1233 +
+ l12*Njq*C1222 +
(l12*Njr+l13*Njq)*C1223 +
((l11+1)*Njr+l13*Njp)*C1213 +
((l11+1)*Njq+l12*Njp)*C1212 +
- (l11+1)*Njp*C1112) +
+ (l11+1)*Njp*C1211) +
(l22+1)*Niq*(l13*Njr*C2233 +
(l12*Njr+l13*Njq)*C2223 +
l12*Njq*C2222 +
((l11+1)*Njr+l13*Njp)*C2213 +
((l11+1)*Njq+l12*Njp)*C2212 +
- (l11+1)*Njp*C1122) +
+ (l11+1)*Njp*C2211) +
l21*Nip*(l13*Njr*C1133 +
(l12*Njr+l13*Njq)*C1123 +
l12*Njq*C1122 +
@@ -919,32 +937,32 @@
((l22+1)*Njr+l23*Njq)*C3323 +
(l21*Njr+l23*Njp)*C3313 +
(l21*Njq+(l22+1)*Njp)*C3312 +
- (l22+1)*Njq*C2233 +
- l21*Njp*C1133) +
- ((l22+1)*Nir+l23*Niq)*(l23*Njr*C3323 +
+ (l22+1)*Njq*C3322 +
+ l21*Njp*C3311) +
+ ((l22+1)*Nir+l23*Niq)*(l23*Njr*C2333 +
((l22+1)*Njr+l23*Njq)*C2323 +
(l21*Njr+l23*Njp)*C2313 +
- (l22+1)*Njq*C2223 +
- (l21*Njq+(l22+1)*Njp)*C1223 +
- l21*Njp*C1123) +
- (l21*Nir+l23*Nip)*(l23*Njr*C3313 +
- ((l22+1)*Njr+l23*Njq)*C2313 +
- (l22+1)*Njq*C2213 +
+ (l22+1)*Njq*C2322 +
+ (l21*Njq+(l22+1)*Njp)*C2312 +
+ l21*Njp*C2311) +
+ (l21*Nir+l23*Nip)*(l23*Njr*C1333 +
+ ((l22+1)*Njr+l23*Njq)*C1323 +
+ (l22+1)*Njq*C1322 +
(l21*Njr+l23*Njp)*C1313 +
- (l21*Njq+(l22+1)*Njp)*C1213 +
- l21*Njp*C1113) +
- (l21*Niq+(l22+1)*Nip)*(l23*Njr*C3312 +
- (l22+1)*Njq*C2212 +
+ (l21*Njq+(l22+1)*Njp)*C1312 +
+ l21*Njp*C1311) +
+ (l21*Niq+(l22+1)*Nip)*(l23*Njr*C1233 +
+ (l22+1)*Njq*C1222 +
((l22+1)*Njr+l23*Njq)*C1223 +
(l21*Njr+l23*Njp)*C1213 +
(l21*Njq+(l22+1)*Njp)*C1212 +
- l21*Njp*C1112) +
+ l21*Njp*C1211) +
(l22+1)*Niq*(l23*Njr*C2233 +
((l22+1)*Njr+l23*Njq)*C2223 +
(l22+1)*Njq*C2222 +
(l21*Njr+l23*Njp)*C2213 +
(l21*Njq+(l22+1)*Njp)*C2212 +
- l21*Njp*C1122) +
+ l21*Njp*C2211) +
l21*Nip*(l23*Njr*C1133 +
((l22+1)*Njr+l23*Njq)*C1123 +
(l22+1)*Njq*C1122 +
@@ -957,69 +975,70 @@
(l32*Njr+(l33+1)*Njq)*C3323 +
(l31*Njr+(l33+1)*Njp)*C3313 +
(l31*Njq+l32*Njp)*C3312 +
- l32*Njq*C2233 +
- l31*Njp*C1133) +
- ((l22+1)*Nir+l23*Niq)*((l33+1)*Njr*C3323 +
+ l32*Njq*C3322 +
+ l31*Njp*C3311) +
+ ((l22+1)*Nir+l23*Niq)*((l33+1)*Njr*C2333 +
(l32*Njr+(l33+1)*Njq)*C2323 +
(l31*Njr+(l33+1)*Njp)*C2313 +
- l32*Njq*C2223 +
- (l31*Njq+l32*Njp)*C1223 +
- l31*Njp*C1123) +
- (l21*Nir+l23*Nip)*((l33+1)*Njr*C3313 +
- (l32*Njr+(l33+1)*Njq)*C2313 +
- l32*Njq*C2213 +
+ l32*Njq*C2322 +
+ (l31*Njq+l32*Njp)*C2312 +
+ l31*Njp*C2311) +
+ (l21*Nir+l23*Nip)*((l33+1)*Njr*C1333 +
+ (l32*Njr+(l33+1)*Njq)*C1323 +
+ l32*Njq*C1322 +
(l31*Njr+(l33+1)*Njp)*C1313 +
- (l31*Njq+l32*Njp)*C1213 +
- l31*Njp*C1113) +
- (l21*Niq+(l22+1)*Nip)*((l33+1)*Njr*C3312 +
- l32*Njq*C2212 +
+ (l31*Njq+l32*Njp)*C1312 +
+ l31*Njp*C1311) +
+ (l21*Niq+(l22+1)*Nip)*((l33+1)*Njr*C1233 +
+ l32*Njq*C1222 +
(l32*Njr+(l33+1)*Njq)*C1223 +
(l31*Njr+(l33+1)*Njp)*C1213 +
(l31*Njq+l32*Njp)*C1212 +
- l31*Njp*C1112) +
+ l31*Njp*C1211) +
(l22+1)*Niq*((l33+1)*Njr*C2233 +
(l32*Njr+(l33+1)*Njq)*C2223 +
l32*Njq*C2222 +
(l31*Njr+(l33+1)*Njp)*C2213 +
(l31*Njq+l32*Njp)*C2212 +
- l31*Njp*C1122) +
+ l31*Njp*C2211) +
l21*Nip*((l33+1)*Njr*C1133 +
(l32*Njr+(l33+1)*Njq)*C1123 +
l32*Njq*C1122 +
(l31*Njr+(l33+1)*Njp)*C1113 +
- (l31*Njq+l32*Njp)*C1112+l31*Njp*C1111);
+ (l31*Njq+l32*Njp)*C1112 +
+ l31*Njp*C1111);
const double Ki2j0 =
(l33+1)*Nir*(l13*Njr*C3333 +
(l12*Njr+l13*Njq)*C3323 +
((l11+1)*Njr+l13*Njp)*C3313 +
((l11+1)*Njq+l12*Njp)*C3312 +
- l12*Njq*C2233 +
- (l11+1)*Njp*C1133) +
- (l32*Nir+(l33+1)*Niq)*(l13*Njr*C3323 +
+ l12*Njq*C3322 +
+ (l11+1)*Njp*C3311) +
+ (l32*Nir+(l33+1)*Niq)*(l13*Njr*C2333 +
(l12*Njr+l13*Njq)*C2323 +
((l11+1)*Njr+l13*Njp)*C2313 +
- l12*Njq*C2223 +
- ((l11+1)*Njq+l12*Njp)*C1223 +
- (l11+1)*Njp*C1123) +
- (l31*Nir+(l33+1)*Nip)*(l13*Njr*C3313 +
- (l12*Njr+l13*Njq)*C2313 +
- l12*Njq*C2213 +
+ l12*Njq*C2322 +
+ ((l11+1)*Njq+l12*Njp)*C2312 +
+ (l11+1)*Njp*C2311) +
+ (l31*Nir+(l33+1)*Nip)*(l13*Njr*C1333 +
+ (l12*Njr+l13*Njq)*C1323 +
+ l12*Njq*C1322 +
((l11+1)*Njr+l13*Njp)*C1313 +
- ((l11+1)*Njq+l12*Njp)*C1213 +
- (l11+1)*Njp*C1113) +
- (l31*Niq+l32*Nip)*(l13*Njr*C3312 +
- l12*Njq*C2212 +
+ ((l11+1)*Njq+l12*Njp)*C1312 +
+ (l11+1)*Njp*C1311) +
+ (l31*Niq+l32*Nip)*(l13*Njr*C1233 +
+ l12*Njq*C1222 +
(l12*Njr+l13*Njq)*C1223 +
((l11+1)*Njr+l13*Njp)*C1213 +
((l11+1)*Njq+l12*Njp)*C1212 +
- (l11+1)*Njp*C1112) +
+ (l11+1)*Njp*C1211) +
l32*Niq*(l13*Njr*C2233 +
(l12*Njr+l13*Njq)*C2223 +
l12*Njq*C2222 +
((l11+1)*Njr+l13*Njp)*C2213 +
((l11+1)*Njq+l12*Njp)*C2212 +
- (l11+1)*Njp*C1122) +
+ (l11+1)*Njp*C2211) +
l31*Nip*(l13*Njr*C1133 +
(l12*Njr+l13*Njq)*C1123 +
l12*Njq*C1122 +
@@ -1032,32 +1051,32 @@
((l22+1)*Njr+l23*Njq)*C3323 +
(l21*Njr+l23*Njp)*C3313 +
(l21*Njq+(l22+1)*Njp)*C3312 +
- (l22+1)*Njq*C2233 +
- l21*Njp*C1133) +
- (l32*Nir+(l33+1)*Niq)*(l23*Njr*C3323 +
+ (l22+1)*Njq*C3322 +
+ l21*Njp*C3311) +
+ (l32*Nir+(l33+1)*Niq)*(l23*Njr*C2333 +
((l22+1)*Njr+l23*Njq)*C2323 +
(l21*Njr+l23*Njp)*C2313 +
- (l22+1)*Njq*C2223 +
- (l21*Njq+(l22+1)*Njp)*C1223 +
- l21*Njp*C1123) +
- (l31*Nir+(l33+1)*Nip)*(l23*Njr*C3313 +
- ((l22+1)*Njr+l23*Njq)*C2313 +
- (l22+1)*Njq*C2213 +
+ (l22+1)*Njq*C2322 +
+ (l21*Njq+(l22+1)*Njp)*C2312 +
+ l21*Njp*C2311) +
+ (l31*Nir+(l33+1)*Nip)*(l23*Njr*C1333 +
+ ((l22+1)*Njr+l23*Njq)*C1323 +
+ (l22+1)*Njq*C1322 +
(l21*Njr+l23*Njp)*C1313 +
- (l21*Njq+(l22+1)*Njp)*C1213 +
- l21*Njp*C1113) +
- (l31*Niq+l32*Nip)*(l23*Njr*C3312 +
- (l22+1)*Njq*C2212 +
+ (l21*Njq+(l22+1)*Njp)*C1312 +
+ l21*Njp*C1311) +
+ (l31*Niq+l32*Nip)*(l23*Njr*C1233 +
+ (l22+1)*Njq*C1222 +
((l22+1)*Njr+l23*Njq)*C1223 +
(l21*Njr+l23*Njp)*C1213 +
(l21*Njq+(l22+1)*Njp)*C1212 +
- l21*Njp*C1112) +
+ l21*Njp*C1211) +
l32*Niq*(l23*Njr*C2233 +
((l22+1)*Njr+l23*Njq)*C2223 +
(l22+1)*Njq*C2222 +
(l21*Njr+l23*Njp)*C2213 +
(l21*Njq+(l22+1)*Njp)*C2212 +
- l21*Njp*C1122) +
+ l21*Njp*C2211) +
l31*Nip*(l23*Njr*C1133 +
((l22+1)*Njr+l23*Njq)*C1123 +
(l22+1)*Njq*C1122 +
@@ -1070,32 +1089,32 @@
(l32*Njr+(l33+1)*Njq)*C3323 +
(l31*Njr+(l33+1)*Njp)*C3313 +
(l31*Njq+l32*Njp)*C3312 +
- l32*Njq*C2233 +
- l31*Njp*C1133) +
- (l32*Nir+(l33+1)*Niq)*((l33+1)*Njr*C3323 +
+ l32*Njq*C3322 +
+ l31*Njp*C3311) +
+ (l32*Nir+(l33+1)*Niq)*((l33+1)*Njr*C2333 +
(l32*Njr+(l33+1)*Njq)*C2323 +
(l31*Njr+(l33+1)*Njp)*C2313 +
- l32*Njq*C2223 +
- (l31*Njq+l32*Njp)*C1223 +
- l31*Njp*C1123) +
- (l31*Nir+(l33+1)*Nip)*((l33+1)*Njr*C3313 +
- (l32*Njr+(l33+1)*Njq)*C2313 +
- l32*Njq*C2213 +
+ l32*Njq*C2322 +
+ (l31*Njq+l32*Njp)*C2312 +
+ l31*Njp*C2311) +
+ (l31*Nir+(l33+1)*Nip)*((l33+1)*Njr*C1333 +
+ (l32*Njr+(l33+1)*Njq)*C1323 +
+ l32*Njq*C1322 +
(l31*Njr+(l33+1)*Njp)*C1313 +
- (l31*Njq+l32*Njp)*C1213 +
- l31*Njp*C1113) +
- (l31*Niq+l32*Nip)*((l33+1)*Njr*C3312 +
- l32*Njq*C2212 +
+ (l31*Njq+l32*Njp)*C1312 +
+ l31*Njp*C1311) +
+ (l31*Niq+l32*Nip)*((l33+1)*Njr*C1233 +
+ l32*Njq*C1222 +
(l32*Njr+(l33+1)*Njq)*C1223 +
(l31*Njr+(l33+1)*Njp)*C1213 +
(l31*Njq+l32*Njp)*C1212 +
- l31*Njp*C1112) +
+ l31*Njp*C1211) +
l32*Niq*((l33+1)*Njr*C2233 +
(l32*Njr+(l33+1)*Njq)*C2223 +
l32*Njq*C2222 +
(l31*Njr+(l33+1)*Njp)*C2213 +
(l31*Njq+l32*Njp)*C2212 +
- l31*Njp*C1122) +
+ l31*Njp*C2211) +
l31*Nip*((l33+1)*Njr*C1133 +
(l32*Njr+(l33+1)*Njq)*C1123 +
l32*Njq*C1122 +
Added: short/3D/PyLith/trunk/libsrc/feassemble/jacobian2d_nonsymm_lgdeform.wxm
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/jacobian2d_nonsymm_lgdeform.wxm (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/jacobian2d_nonsymm_lgdeform.wxm 2010-03-17 01:03:36 UTC (rev 16429)
@@ -0,0 +1,53 @@
+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 0.8.2 ] */
+
+/* [wxMaxima: input start ] */
+Bnl: matrix([N11,0,N21,0,N31,0],
+[N12,0,N22,0,N32,0],
+[0,N11,0,N21,0,N31],
+[0,N12,0,N22,0,N32]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+S: matrix([S11,S12,0,0],
+[S12,S22,0,0],
+[0,0,S11,S12],
+[0,0,S12,S22]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+x: transpose(Bnl) . S . Bnl;
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Bi: matrix([Nip*(1+l11), l21*Nip],
+[l12*Niq,(1+l22)*Niq],
+[(1+l11)*Niq+l12*Nip,(1+l22)*Nip+l21*Niq]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Bj: matrix([Njp*(1+l11), l21*Njp],
+[l12*Njq,(1+l22)*Njq],
+[(1+l11)*Njq+l12*Njp,(1+l22)*Njp+l21*Njq]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+C: matrix([C1111, C1122, C1112],
+[C2211, C2222, C2212],
+[C1211, C1222, C1212]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+y: transpose(Bi) . C . Bj;
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Svec: matrix([s11, s22, s12]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+r: transpose(Bi) . transpose(Svec);
+/* [wxMaxima: input end ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$
Added: short/3D/PyLith/trunk/libsrc/feassemble/jacobian3d_nonsymm_lgdeform.wxm
===================================================================
--- short/3D/PyLith/trunk/libsrc/feassemble/jacobian3d_nonsymm_lgdeform.wxm (rev 0)
+++ short/3D/PyLith/trunk/libsrc/feassemble/jacobian3d_nonsymm_lgdeform.wxm 2010-03-17 01:03:36 UTC (rev 16429)
@@ -0,0 +1,84 @@
+/* [wxMaxima batch file version 1] [ DO NOT EDIT BY HAND! ]*/
+/* [ Created with wxMaxima version 0.8.2 ] */
+
+/* [wxMaxima: input start ] */
+Bnli: matrix([Nip,0,0],
+[Niq,0,0],
+[Nir,0,0],
+[0,Nip,0],
+[0,Niq,0],
+[0,Nir,0],
+[0,0,Nip],
+[0,0,Niq],
+[0,0,Nir]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Bnlj: matrix([Njp,0,0],
+[Njq,0,0],
+[Njr,0,0],
+[0,Njp,0],
+[0,Njq,0],
+[0,Njr,0],
+[0,0,Njp],
+[0,0,Njq],
+[0,0,Njr]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+S: matrix([s11,s12,s13,0,0,0,0,0,0],
+[s12,s22,s23,0,0,0,0,0,0],
+[s13,s23,s33,0,0,0,0,0,0],
+[0,0,0,s11,s12,s13,0,0,0],
+[0,0,0,s12,s22,s23,0,0,0],
+[0,0,0,s13,s23,s33,0,0,0],
+[0,0,0,0,0,0,s11,s12,s13],
+[0,0,0,0,0,0,s12,s22,s23],
+[0,0,0,0,0,0,s13,s23,s33]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+x: transpose(Bnli) . S . Bnlj;
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Bi: matrix([Nip*(1+l11), l21*Nip, l31*Nip],
+[l12*Niq,(1+l22)*Niq, l32*Niq],
+[l13*Nir, l23*Nir, (1+l33)*Nir],
+[(1+l11)*Niq+l12*Nip,(1+l22)*Nip+l21*Niq,l31*Niq+l32*Nip],
+[l12*Nir+l13*Niq, (1+l22)*Nir+l23*Niq, (1+l33)*Niq+l32*Nir],
+[(1+l11)*Nir+l13*Nip, l21*Nir+l23*Nip, (1+l33)*Nip+l31*Nir]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Bj: matrix([Njp*(1+l11), l21*Njp, l31*Njp],
+[l12*Njq,(1+l22)*Njq, l32*Njq],
+[l13*Njr, l23*Njr, (1+l33)*Njr],
+[(1+l11)*Njq+l12*Njp,(1+l22)*Njp+l21*Njq,l31*Njq+l32*Njp],
+[l12*Njr+l13*Njq, (1+l22)*Njr+l23*Njq, (1+l33)*Njq+l32*Njr],
+[(1+l11)*Njr+l13*Njp, l21*Njr+l23*Njp, (1+l33)*Njp+l31*Njr]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+C: matrix([C1111, C1122, C1133, C1112, C1123, C1113],
+[C2211, C2222, C2233, C2212, C2223, C2213],
+[C3311, C3322, C3333, C3312, C3323, C3313],
+[C1211, C1222, C1233, C1212, C1223, C1213],
+[C2311, C2322, C2333, C2312, C2323, C2313],
+[C1311, C1322, C1333, C1312, C1323, C1313]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+y: transpose(Bi) . C . Bj;
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+Svec: matrix([s11,s22,s33,s12,s23,s13]);
+/* [wxMaxima: input end ] */
+
+/* [wxMaxima: input start ] */
+r: transpose(Bi) . transpose(Svec);
+/* [wxMaxima: input end ] */
+
+/* Maxima can't load/batch files which end with a comment! */
+"Created with wxMaxima"$
Modified: short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.hh
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.hh 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.hh 2010-03-17 01:03:36 UTC (rev 16429)
@@ -446,7 +446,21 @@
double _scalarProduct(const double* tensor1,
const double* tensor2) const;
+ /** Compute tensor mean, assuming vector form of a tensor.
+ *
+ * @param vec Tensor represented as a vector.
+ */
+ double _calcMean(const double* vec) const;
+ /** Compute deviatoric components, assuming vector form of a tensor.
+ *
+ * @param vec Tensor represented as a vector.
+ * @param vecMean Tensor mean.
+ */
+ double _calcDeviatoric(const double* vec,
+ const double vecMean) const;
+
+
// PRIVATE MEMBERS ////////////////////////////////////////////////////
private :
Modified: short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.icc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.icc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/DruckerPrager3D.icc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -27,6 +27,28 @@
_dt = dt;
} // timeStep
+// Compute mean stress/strain from vector.
+inline
+double
+pylith::materials::DruckerPragerEP3D::_calcMean(const double* vec) {
+ const double vecMean = (vec[0] + vec[1] + vec[2])/3.0;
+ return vecMean;
+} // _calcMean
+
+// Compute deviatoric stress/strain from vector and mean value.
+inline
+double
+pylith::materials::DruckerPragerEP3D::_calcDeviatoric(const double* vec,
+ const double vecMean) {
+ const double deviatoric[] = {vec[0] - vecMean,
+ vec[1] - vecMean,
+ vec[2] - vecMean,
+ vec[3],
+ vec[4],
+ vec[5]};
+ return deviatoric;
+} // _calcDeviatoric
+
// Compute stress tensor from parameters.
inline
void
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticIsotropic3D.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -39,7 +39,7 @@
const int tensorSize = 6;
// Number of elastic constants (for general 3-D elastic material)
- const int numElasticConsts = 21;
+ const int numElasticConsts = 36;
// Number of physical properties.
const int numProperties = 3;
@@ -284,21 +284,36 @@
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = mu2; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = mu2; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = mu2; // C1313
+ elasticConsts[ 6] = lambda; // C2211
+ elasticConsts[ 7] = lambda2mu; // C2222
+ elasticConsts[ 8] = lambda; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = lambda; // C3311
+ elasticConsts[13] = lambda; // C3322
+ elasticConsts[14] = lambda2mu; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = mu2; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = mu2; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = mu2; // C1313
PetscLogFlops(2);
} // _calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStrain.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -39,7 +39,7 @@
const int tensorSize = 3;
// Number of elastic constants (for general 3-D elastic material)
- const int numElasticConsts = 6;
+ const int numElasticConsts = 9;
// Number of physical properties.
const int numProperties = 3;
@@ -292,9 +292,12 @@
elasticConsts[0] = lambda2mu; // C1111
elasticConsts[1] = lambda; // C1122
elasticConsts[2] = 0; // C1112
- elasticConsts[3] = lambda2mu; // C2222
- elasticConsts[4] = 0; // C2212
- elasticConsts[5] = mu2; // C1212
+ elasticConsts[3] = lambda; // C2211
+ elasticConsts[4] = lambda2mu; // C2222
+ elasticConsts[5] = 0; // C2212
+ elasticConsts[6] = 0; // C1211
+ elasticConsts[7] = 0; // C1222
+ elasticConsts[8] = mu2; // C1212
PetscLogFlops(2);
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/ElasticPlaneStress.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -38,8 +38,8 @@
// Number of entries in stress tensor.
const int tensorSize = 3;
- // Number of elastic constants (for general 3-D elastic material)
- const int numElasticConsts = 6;
+ // Number of elastic constants (for general 2-D elastic material)
+ const int numElasticConsts = 9;
// Number of physical properties.
const int numProperties = 3;
@@ -295,9 +295,12 @@
elasticConsts[0] = c11; // C1111
elasticConsts[1] = mu2 * lambda / lambda2mu; // C1122
elasticConsts[2] = 0; // C1112
- elasticConsts[3] = c11; // C2222
- elasticConsts[4] = 0; // C2212
- elasticConsts[5] = mu2; // C1212
+ elasticConsts[3] = elasticConsts[1]; // C2211
+ elasticConsts[4] = c11; // C2222
+ elasticConsts[5] = 0; // C2212
+ elasticConsts[6] = 0; // C1211
+ elasticConsts[7] = 0; // C1222
+ elasticConsts[8] = mu2; // C1212
PetscLogFlops(8);
} // calcElasticConsts
Modified: short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/GenMaxwellIsotropic3D.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -44,7 +44,7 @@
const int tensorSize = 6;
/// Number of entries in derivative of elasticity matrix.
- const int numElasticConsts = 21;
+ const int numElasticConsts = 36;
/// Number of physical properties.
const int numProperties = 5;
@@ -596,23 +596,38 @@
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = mu2; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = mu2; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = mu2; // C1313
+ elasticConsts[ 6] = lambda; // C2211
+ elasticConsts[ 7] = lambda2mu; // C2222
+ elasticConsts[ 8] = lambda; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = lambda; // C3311
+ elasticConsts[13] = lambda; // C3322
+ elasticConsts[14] = lambda2mu; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = mu2; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = mu2; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = mu2; // C1313
- PetscLogFlops(4);
+ PetscLogFlops(2);
} // _calcElasticConstsElastic
// ----------------------------------------------------------------------
@@ -670,27 +685,42 @@
double elasFrac = 1.0 - visFrac;
double shearFac = mu*(elasFrac + visFac)/3.0;
- elasticConsts[ 0] = bulkModulus + 4.0*shearFac; // C1111
- elasticConsts[ 1] = bulkModulus - 2.0*shearFac; // C1122
+ elasticConsts[ 0] = bulkModulus + 4.0 * shearFac; // C1111
+ elasticConsts[ 1] = bulkModulus - 2.0 * shearFac; // C1122
elasticConsts[ 2] = elasticConsts[1]; // C1133
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = elasticConsts[0]; // C2222
- elasticConsts[ 7] = elasticConsts[1]; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = elasticConsts[0]; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = 6.0 * shearFac; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = elasticConsts[15]; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = elasticConsts[15]; // C1313
+ elasticConsts[ 6] = elasticConsts[1]; // C2211
+ elasticConsts[ 7] = elasticConsts[0]; // C2222
+ elasticConsts[ 8] = elasticConsts[1]; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = elasticConsts[1]; // C3311
+ elasticConsts[13] = elasticConsts[1]; // C3322
+ elasticConsts[14] = elasticConsts[0]; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = 6.0 * shearFac; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = elasticConsts[21]; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = elasticConsts[21]; // C1313
#if 0 // DEBUGGING
std::cout << "_calcElasticConstsViscoelastic" << std::endl;
@@ -775,7 +805,7 @@
const int initialStressSize,
const double* initialStrain,
const int initialStrainSize)
-{ // _updateStateVarsElastic
+{ // _updateStateVarsViscoelastic
assert(0 != stateVars);
assert(_numVarsQuadPt == numStateVars);
assert(0 != properties);
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellIsotropic3D.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -40,7 +40,7 @@
const int tensorSize = 6;
/// Number of entries in derivative of elasticity matrix.
- const int numElasticConsts = 21;
+ const int numElasticConsts = 36;
/// Number of physical properties.
const int numProperties = 4;
@@ -503,23 +503,38 @@
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = mu2; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = mu2; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = mu2; // C1313
+ elasticConsts[ 6] = lambda; // C2211
+ elasticConsts[ 7] = lambda2mu; // C2222
+ elasticConsts[ 8] = lambda; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = lambda; // C3311
+ elasticConsts[13] = lambda; // C3322
+ elasticConsts[14] = lambda2mu; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = mu2; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = mu2; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = mu2; // C1313
- PetscLogFlops(4);
+ PetscLogFlops(2);
} // _calcElasticConstsElastic
// ----------------------------------------------------------------------
@@ -563,27 +578,43 @@
double dq = ViscoelasticMaxwell::viscousStrainParam(_dt, maxwellTime);
const double visFac = mu * dq / 3.0;
+
elasticConsts[ 0] = bulkModulus + 4.0 * visFac; // C1111
elasticConsts[ 1] = bulkModulus - 2.0 * visFac; // C1122
elasticConsts[ 2] = elasticConsts[1]; // C1133
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = elasticConsts[0]; // C2222
- elasticConsts[ 7] = elasticConsts[1]; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = elasticConsts[0]; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = 6.0 * visFac; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = elasticConsts[15]; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = elasticConsts[15]; // C1313
+ elasticConsts[ 6] = elasticConsts[1]; // C2211
+ elasticConsts[ 7] = elasticConsts[0]; // C2222
+ elasticConsts[ 8] = elasticConsts[1]; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = elasticConsts[1]; // C3311
+ elasticConsts[13] = elasticConsts[1]; // C3322
+ elasticConsts[14] = elasticConsts[0]; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = 6.0 * visFac; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = elasticConsts[21]; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = elasticConsts[21]; // C1313
PetscLogFlops(10);
} // _calcElasticConstsViscoelastic
Modified: short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/MaxwellPlaneStrain.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -40,7 +40,7 @@
const int tensorSize = 3;
/// Number of entries in derivative of elasticity matrix.
- const int numElasticConsts = 6;
+ const int numElasticConsts = 9;
/// Number of physical properties.
const int numProperties = 4;
@@ -477,9 +477,12 @@
elasticConsts[ 0] = lambda2mu; // C1111
elasticConsts[ 1] = lambda; // C1122
elasticConsts[ 2] = 0; // C1112
- elasticConsts[ 3] = lambda2mu; // C2222
- elasticConsts[ 4] = 0; // C2212
- elasticConsts[ 5] = mu2; // C1212
+ elasticConsts[ 3] = lambda; // C2211
+ elasticConsts[ 4] = lambda2mu; // C2222
+ elasticConsts[ 5] = 0; // C2212
+ elasticConsts[ 6] = 0; // C1211
+ elasticConsts[ 7] = 0; // C1222
+ elasticConsts[ 8] = mu2; // C1212
PetscLogFlops(2);
} // _calcElasticConstsElastic
@@ -528,9 +531,12 @@
elasticConsts[ 0] = bulkModulus + 4.0 * visFac; // C1111
elasticConsts[ 1] = bulkModulus - 2.0 * visFac; // C1122
elasticConsts[ 2] = 0; // C1112
- elasticConsts[ 3] = elasticConsts[0]; // C2222
- elasticConsts[ 4] = 0; // C2212
- elasticConsts[ 5] = 6.0 * visFac; // C1212
+ elasticConsts[ 3] = elasticConsts[1]; // C2211
+ elasticConsts[ 4] = elasticConsts[0]; // C2222
+ elasticConsts[ 5] = 0; // C2212
+ elasticConsts[ 6] = 0; // C1211
+ elasticConsts[ 7] = 0; // C1222
+ elasticConsts[ 8] = 6.0 * visFac; // C1212
PetscLogFlops(10);
} // _calcElasticConstsViscoelastic
Modified: short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc
===================================================================
--- short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2010-03-16 00:28:24 UTC (rev 16428)
+++ short/3D/PyLith/trunk/libsrc/materials/PowerLaw3D.cc 2010-03-17 01:03:36 UTC (rev 16429)
@@ -42,7 +42,7 @@
const int tensorSize = 6;
/// Number of entries in derivative of elasticity matrix.
- const int numElasticConsts = 21;
+ const int numElasticConsts = 36;
/// Number of physical properties.
const int numProperties = 6;
@@ -807,21 +807,36 @@
elasticConsts[ 3] = 0; // C1112
elasticConsts[ 4] = 0; // C1123
elasticConsts[ 5] = 0; // C1113
- elasticConsts[ 6] = lambda2mu; // C2222
- elasticConsts[ 7] = lambda; // C2233
- elasticConsts[ 8] = 0; // C2212
- elasticConsts[ 9] = 0; // C2223
- elasticConsts[10] = 0; // C2213
- elasticConsts[11] = lambda2mu; // C3333
- elasticConsts[12] = 0; // C3312
- elasticConsts[13] = 0; // C3323
- elasticConsts[14] = 0; // C3313
- elasticConsts[15] = mu2; // C1212
- elasticConsts[16] = 0; // C1223
- elasticConsts[17] = 0; // C1213
- elasticConsts[18] = mu2; // C2323
- elasticConsts[19] = 0; // C2313
- elasticConsts[20] = mu2; // C1313
+ elasticConsts[ 6] = lambda; // C2211
+ elasticConsts[ 7] = lambda2mu; // C2222
+ elasticConsts[ 8] = lambda; // C2233
+ elasticConsts[ 9] = 0; // C2212
+ elasticConsts[10] = 0; // C2223
+ elasticConsts[11] = 0; // C2213
+ elasticConsts[12] = lambda; // C3311
+ elasticConsts[13] = lambda; // C3322
+ elasticConsts[14] = lambda2mu; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = mu2; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = mu2; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = mu2; // C1313
PetscLogFlops(2);
} // _calcElasticConstsElastic
@@ -1041,28 +1056,42 @@
/// Compute tangent matrix.
// Form elastic constants.
- elasticConsts[ 0] = bulkModulus + 2.0 * dStress11dStrain11/3.0; // C1111
- elasticConsts[ 1] = bulkModulus - dStress11dStrain11/3.0; // C1122
+ elasticConsts[ 0] = bulkModulus + 2.0 * dStress11dStrain11/3.0; // C1111
+ elasticConsts[ 1] = bulkModulus - dStress11dStrain11/3.0; // C1122
elasticConsts[ 2] = elasticConsts[ 1]; // C1133
- elasticConsts[ 3] = 0.0; // C1112
- elasticConsts[ 4] = 0.0; // C1123
- elasticConsts[ 5] = 0.0; // C1113
- elasticConsts[ 6] = bulkModulus + 2.0 * dStress22dStrain22/3.0; // C2222
- elasticConsts[ 7] = bulkModulus - dStress22dStrain22/3.0; // C2233
- elasticConsts[ 8] = 0.0; // C2212
- elasticConsts[ 9] = 0.0; // C2223
- elasticConsts[10] = 0.0; // C2213
- elasticConsts[11] = bulkModulus + 2.0 * dStress33dStrain33/3.0; // C3333
- elasticConsts[12] = 0.0; // C3312
- elasticConsts[13] = 0.0; // C3323
- elasticConsts[14] = 0.0; // C3313
- elasticConsts[15] = dStress12dStrain12; // C1212
- elasticConsts[16] = 0.0; // C1223
- elasticConsts[17] = 0.0; // C1213
- elasticConsts[18] = dStress23dStrain23; // C2323
- elasticConsts[19] = 0.0; // C2313
- elasticConsts[20] = dStress13dStrain13; // C1313
-
+ elasticConsts[ 3] = 0.0; // C1112
+ elasticConsts[ 4] = 0.0; // C1123
+ elasticConsts[ 5] = 0.0; // C1113
+ elasticConsts[ 6] = elasticConsts[ 1]; // C2211
+ elasticConsts[ 7] = bulkModulus + 2.0 * dStress22dStrain22/3.0; // C2222
+ elasticConsts[ 8] = bulkModulus - dStress22dStrain22/3.0; // C2233
+ elasticConsts[ 9] = 0.0; // C2212
+ elasticConsts[10] = 0.0; // C2223
+ elasticConsts[11] = 0.0; // C2213
+ elasticConsts[12] = elasticConsts[ 1]; // C3311
+ elasticConsts[13] = elasticConsts[ 8]; // C3322
+ elasticConsts[14] = bulkModulus + 2.0 * dStress33dStrain33/3.0; // C3333
+ elasticConsts[15] = 0; // C3312
+ elasticConsts[16] = 0; // C3323
+ elasticConsts[17] = 0; // C3313
+ elasticConsts[18] = 0; // C1211
+ elasticConsts[19] = 0; // C1222
+ elasticConsts[20] = 0; // C1233
+ elasticConsts[21] = dStress12dStrain12; // C1212
+ elasticConsts[22] = 0; // C1223
+ elasticConsts[23] = 0; // C1213
+ elasticConsts[24] = 0; // C2311
+ elasticConsts[25] = 0; // C2322
+ elasticConsts[26] = 0; // C2333
+ elasticConsts[27] = 0; // C2312
+ elasticConsts[28] = dStress23dStrain23; // C2323
+ elasticConsts[29] = 0; // C2313
+ elasticConsts[30] = 0; // C1311
+ elasticConsts[31] = 0; // C1322
+ elasticConsts[32] = 0; // C1333
+ elasticConsts[33] = 0; // C1312
+ elasticConsts[34] = 0; // C1323
+ elasticConsts[35] = dStress13dStrain13; // C1313
PetscLogFlops(114);
} // else
} // _calcElasticConstsViscoelastic
More information about the CIG-COMMITS
mailing list