use of maspack.matrix.SymmetricMatrix3d in project artisynth_core by artisynth.
the class FungMaterial method computeTangent.
public void computeTangent(Matrix6d c, SymmetricMatrix3d stress, SolidDeformation def, Matrix3d Q, FemMaterial baseMat) {
c.setZero();
double[] K = new double[3];
double[] L = new double[3];
Vector3d[] a0 = { new Vector3d(), new Vector3d(), new Vector3d() };
Vector3d[] a = { new Vector3d(), new Vector3d(), new Vector3d() };
Vector3d vtmp = new Vector3d();
// Matrix3d Q = Matrix3d.IDENTITY;
// if (dt.getFrame() != null) {
// Q = dt.getFrame();
// }
SymmetricMatrix3d[] A = { new SymmetricMatrix3d(), new SymmetricMatrix3d(), new SymmetricMatrix3d() };
Matrix3d tmpMatrix = new Matrix3d();
Matrix3d tmpMatrix2 = new Matrix3d();
SymmetricMatrix3d tmpSymmMatrix = new SymmetricMatrix3d();
Matrix6d tmpMatrix6d = new Matrix6d();
Matrix6d cFung = new Matrix6d();
// Evaluate Lame coefficients
mu[0] = myMU1;
mu[1] = myMU2;
mu[2] = myMU3;
lam[0][0] = myL11;
lam[0][1] = myL12;
lam[0][2] = myL31;
lam[1][0] = myL12;
lam[1][1] = myL22;
lam[1][2] = myL23;
lam[2][0] = myL31;
lam[2][1] = myL23;
lam[2][2] = myL33;
double J = def.getDetF();
double avgp = def.getAveragePressure();
// Calculate deviatoric left Cauchy-Green tensor
def.computeDevLeftCauchyGreen(myB);
// Calculate deviatoric right Cauchy-Green tensor
def.computeDevRightCauchyGreen(myC);
// calculate square of C
myC2.mulTransposeLeft(myC);
Matrix3d mydevF = new Matrix3d(def.getF());
mydevF.scale(Math.pow(J, -1.0 / 3.0));
for (int i = 0; i < 3; i++) {
// Copy the texture direction in the reference configuration to a0
a0[i].x = Q.get(0, i);
a0[i].y = Q.get(1, i);
a0[i].z = Q.get(2, i);
vtmp.mul(myC, a0[i]);
K[i] = a0[i].dot(vtmp);
vtmp.mul(myC2, a0[i]);
L[i] = a0[i].dot(vtmp);
a[i].mul(mydevF, a0[i]);
a[i].scale(Math.pow(K[i], -0.5));
A[i].set(a[i].x * a[i].x, a[i].y * a[i].y, a[i].z * a[i].z, a[i].x * a[i].y, a[i].x * a[i].z, a[i].y * a[i].z);
}
// Evaluate exp(Q)
double eQ = 0.0;
for (int i = 0; i < 3; i++) {
eQ += 2.0 * mu[i] * (L[i] - 2.0 * K[i] + 1.0);
for (int j = 0; j < 3; j++) eQ += lam[i][j] * (K[i] - 1.0) * (K[j] - 1.0);
}
eQ = Math.exp(eQ / (4. * myCC));
// Evaluate the distortional part of the Cauchy stress
SymmetricMatrix3d sd = new SymmetricMatrix3d();
SymmetricMatrix3d bmi = new SymmetricMatrix3d();
bmi.sub(myB, SymmetricMatrix3d.IDENTITY);
for (int i = 0; i < 3; i++) {
// sd += mu[i]*K[i]*(A[i]*bmi + bmi*A[i]);
tmpMatrix.mul(A[i], bmi);
tmpMatrix2.mul(bmi, A[i]);
tmpMatrix.add(tmpMatrix2);
tmpMatrix.scale(mu[i] * K[i]);
tmpSymmMatrix.setSymmetric(tmpMatrix);
sd.add(tmpSymmMatrix);
for (int j = 0; j < 3; j++) {
// s += lam[i][j]*((K[i]-1)*K[j]*A[j]+(K[j]-1)*K[i]*A[i])/2.;
sd.scaledAdd(lam[i][j] / 2.0 * (K[i] - 1.0) * K[j], A[j]);
sd.scaledAdd(lam[i][j] / 2.0 * (K[j] - 1.0) * K[i], A[i]);
}
addTensorProduct4(cFung, mu[i] * K[i], A[i], myB);
// C += mu[i]*K[i]*dyad4s(A[i],b);
for (int j = 0; j < 3; j++) {
TensorUtils.addSymmetricTensorProduct(cFung, lam[i][j] * K[i] * K[j] / 2.0, A[i], A[j]);
// C += lam[i][j]*K[i]*K[j]*dyad1s(A[i],A[j])/2.;
}
}
sd.scale(eQ / (2.0 * J));
// This is the distortional part of the elasticity tensor
// C = (eQ / J)*C + (2.0*J/myc/eQ)*dyad1s(sd);
cFung.scale(eQ / J);
// Factor of two diff between FEBio and Artisynth tensor utility Taken into account with scalefactor
TensorUtils.addSymmetricTensorProduct(cFung, J / myCC / eQ, sd, sd);
// This is the final value of the elasticity tensor
// tens4ds IxI = dyad1s(I);
// tens4ds I4 = dyad4s(I);
double cTrace = cFung.m00 + cFung.m11 + cFung.m22 + cFung.m33 + cFung.m44 + cFung.m55;
// C += - 1./3.*(ddots(C,IxI) - IxI*(C.tr()/3.))
// + 2./3.*((I4-IxI/3.)*sd.tr()-dyad1s(sd.dev(),I));
TensorUtils.addScaledIdentityProduct(tmpMatrix6d, -1.0 / 3.0);
cFung.add(ddots(cFung, tmpMatrix6d));
TensorUtils.addScaledIdentityProduct(cFung, 1.0 / 9.0 * cTrace);
TensorUtils.addTensorProduct4(cFung, 2.0 / 3.0 * sd.trace(), // check
Matrix3d.IDENTITY);
TensorUtils.addScaledIdentityProduct(cFung, -2.0 / 9.0 * sd.trace());
sd.deviator();
tmpMatrix6d.setZero();
TensorUtils.addSymmetricIdentityProduct(tmpMatrix6d, sd);
tmpMatrix6d.scale(-2.0 / 3.0);
cFung.add(tmpMatrix6d);
tmpMatrix6d.setZero();
cFung.setLowerToUpper();
c.add(cFung);
c.m00 += -avgp;
c.m11 += -avgp;
c.m22 += -avgp;
c.m01 += avgp;
c.m02 += avgp;
c.m12 += avgp;
c.m33 += -avgp;
c.m44 += -avgp;
c.m55 += -avgp;
c.setLowerToUpper();
}
use of maspack.matrix.SymmetricMatrix3d in project artisynth_core by artisynth.
the class IncompNeoHookeanMaterial method clone.
public IncompNeoHookeanMaterial clone() {
IncompNeoHookeanMaterial mat = (IncompNeoHookeanMaterial) super.clone();
mat.myB = new SymmetricMatrix3d();
return mat;
}
use of maspack.matrix.SymmetricMatrix3d in project artisynth_core by artisynth.
the class IncompNeoHookeanMaterial method main.
public static void main(String[] args) {
IncompNeoHookeanMaterial mat = new IncompNeoHookeanMaterial();
SolidDeformation def = new SolidDeformation();
Matrix3d Q = new Matrix3d();
def.setF(new Matrix3d(1, 3, 5, 2, 1, 4, 6, 1, 2));
Matrix6d D = new Matrix6d();
SymmetricMatrix3d sig = new SymmetricMatrix3d();
mat.setShearModulus(10);
mat.computeStress(sig, def, Q, null);
mat.computeTangent(D, sig, def, Q, null);
System.out.println("sig=\n" + sig.toString("%12.6f"));
System.out.println("D=\n" + D.toString("%12.6f"));
}
use of maspack.matrix.SymmetricMatrix3d in project artisynth_core by artisynth.
the class BlemkerMuscle method main.
public static void main(String[] args) {
BlemkerMuscle mat = new BlemkerMuscle();
SolidDeformation def = new SolidDeformation();
def.setF(new Matrix3d(1, 3, 5, 2, 1, 4, 6, 1, 2));
Matrix6d D = new Matrix6d();
SymmetricMatrix3d sig = new SymmetricMatrix3d();
FemMaterial baseMat = new MooneyRivlinMaterial();
Vector3d a = new Vector3d(1, 0, 0);
// a.setRandom();
mat.computeStress(sig, 1.0, a, def, baseMat);
// def.setStress (sig);
mat.computeTangent(D, sig, 1.0, a, def, baseMat);
System.out.println("sig=\n" + sig.toString("%12.6f"));
System.out.println("D=\n" + D.toString("%12.6f"));
}
use of maspack.matrix.SymmetricMatrix3d in project artisynth_core by artisynth.
the class OgdenMaterial method clone.
public OgdenMaterial clone() {
OgdenMaterial mat = (OgdenMaterial) super.clone();
mat.myB = new SymmetricMatrix3d();
return mat;
}
Aggregations