Examples with Matrix3d maspack.matrix.Matrix3d used on opensource projects

Example 16 with Matrix3d

use of maspack.matrix.Matrix3d in project artisynth_core by artisynth.

the class OgdenMaterial method computeTangent.

// See JC Simo and RL Taylor, Quasi-Incompressible Finite Elasticity in
// Principal Stretches.  Continuum Basis and Numerical Algorithms, Comp
// Methods Appl Mech Eng, 85, 273-310, 1991
public void computeTangent(Matrix6d c, SymmetricMatrix3d stress, SolidDeformation def, Matrix3d Q, FemMaterial baseMat) {
    double J = def.getDetF();
    double[][] principalStretchDevPow;
    principalStretchDevPow = new double[3][6];
    // Calculate left Cauchy-Green tensor
    def.computeLeftCauchyGreen(myB);
    // Calculate square of B
    myB2.mulTransposeLeft(myB);
    // average element pressure
    double p = def.getAveragePressure();
    Vector3d principalStretch = new Vector3d();
    Vector3d principalStretch2 = new Vector3d();
    Vector3d principalStretchDev = new Vector3d();
    Matrix3d principalDirection = new Matrix3d();
    // Calculate principal stretches and principal directions
    myB.getEigenValues(principalStretch2, principalDirection);
    for (int i = 0; i < 3; i++) {
        principalStretch.set(i, Math.sqrt(principalStretch2.get(i)));
    }
    // Calculate deviatoric principal stretches
    principalStretchDev.scale(Math.pow(J, -1.0 / 3.0), principalStretch);
    // Calculate 1st and 3rd strain invariants
    double I1 = (principalStretch2.get(0) + principalStretch2.get(1) + principalStretch2.get(2));
    double I3 = (principalStretch2.get(0) * principalStretch2.get(1) * principalStretch2.get(2));
    // Raise deviatoric stretches to power of alpha1, alpha2, etc
    for (int i = 0; i < 3; i++) {
        for (int n = 0; n < Nmax; n++) {
            if (myMu[n] != 0) {
                principalStretchDevPow[i][n] = Math.pow(principalStretchDev.get(i), myAlpha[n]);
            }
        }
    }
    SymmetricMatrix3d ma = new SymmetricMatrix3d();
    SymmetricMatrix3d mb = new SymmetricMatrix3d();
    // Constant for determining differences in principal stretches
    double smallNum = 1e-8;
    // Initialise elasticity tensor
    c.setZero();
    // Consider the three distinct cases
    if ((Math.abs(principalStretch2.get(0) - principalStretch2.get(1)) < smallNum) && (Math.abs(principalStretch2.get(0) - principalStretch2.get(2)) < smallNum)) {
        // All three principal stretches are the same
        // gamma_AB in Eq. 2.66 of Simo and Taylor (1991)
        double g = 0.0;
        for (int n = 0; n < Nmax; n++) {
            if (myMu[n] != 0) {
                g += myMu[n] * principalStretchDevPow[0][n];
            }
        }
        // Eq. 2.71 of Simo and Taylor (1991)
        TensorUtils.addScaledIdentity(c, g / J);
        TensorUtils.addScaledIdentityProduct(c, -g / J / 3.0);
    } else if ((Math.abs(principalStretch2.get(0) - principalStretch2.get(1)) >= smallNum) && (Math.abs(principalStretch2.get(0) - principalStretch2.get(2)) >= smallNum) && (Math.abs(principalStretch2.get(1) - principalStretch2.get(2)) >= smallNum)) {
        for (int i = 0; i < 3; i++) {
            int j = (i + 1) % 3;
            int k = (j + 1) % 3;
            // D prime i - Eq. 2.46a of Simo and Taylor (1991)
            double Dpi = 8.0 * Math.pow(principalStretch.get(i), 3.0) - 2.0 * I1 * principalStretch.get(i) - 2.0 * I3 * Math.pow(principalStretch.get(i), -3.0);
            // D i - Eq. 2.16a of Simo and Taylor (1991)
            double Di = (principalStretch2.get(i) - principalStretch2.get(j)) * (principalStretch2.get(i) - principalStretch2.get(k));
            // the matrix mi - Eq. 2.15 of Simo and Taylor (1991)
            ma.set(myB2);
            ma.scaledAdd(-principalStretch2.get(k), myB);
            ma.scaledAdd(-principalStretch2.get(j), myB);
            ma.scaledAdd(principalStretch2.get(j) * principalStretch2.get(k), SymmetricMatrix3d.IDENTITY);
            ma.scale(1.0 / Di);
            // beta_i in Eq. 2.65 of Simo and Taylor (1991)
            double beta = 0.0;
            for (int n = 0; n < Nmax; n++) {
                if (myMu[n] != 0) {
                    beta += myMu[n] / myAlpha[n] * (principalStretchDevPow[i][n] - (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 3.0);
                }
            }
            // Calculate dgm term in Eq 2.68 of Simo and Taylor (1991)
            TensorUtils.addTensorProduct4(c, 2.0 * beta / J / Di, myB);
            TensorUtils.addTensorProduct(c, -2.0 * beta / J / Di, myB);
            TensorUtils.addScaledIdentityProduct(c, I3 * 2.0 * beta / J / Di / principalStretch2.get(i));
            TensorUtils.addScaledIdentity(c, -I3 * 2.0 * beta / J / Di / principalStretch2.get(i));
            TensorUtils.addSymmetricTensorProduct(c, 2.0 * beta / J / Di * principalStretch2.get(i), myB, ma);
            TensorUtils.addTensorProduct(c, -1.0 * beta / J / Di * Dpi * principalStretch.get(i), ma);
            TensorUtils.addSymmetricTensorProduct(c, -2.0 * beta / J / Di * I3 / principalStretch2.get(i), SymmetricMatrix3d.IDENTITY, ma);
            for (int n = 0; n < 3; n++) {
                j = (n + 1) % 3;
                k = (j + 1) % 3;
                Di = (principalStretch2.get(n) - principalStretch2.get(j)) * (principalStretch2.get(n) - principalStretch2.get(k));
                // the matrix mi - Eq. 2.15 of Simo and Taylor (1991)
                mb.set(myB2);
                mb.scaledAdd(-principalStretch2.get(k), myB);
                mb.scaledAdd(-principalStretch2.get(j), myB);
                mb.scaledAdd(principalStretch2.get(j) * principalStretch2.get(k), SymmetricMatrix3d.IDENTITY);
                mb.scale(1.0 / Di);
                // gamma_AB in Eq. 2.66 of Simo and Taylor (1991)
                double gab = 0.0;
                if (i == n) {
                    for (int m = 0; m < Nmax; m++) {
                        if (myMu[m] != 0) {
                            gab += myMu[m] * (principalStretchDevPow[i][m] / 3.0 + (principalStretchDevPow[0][m] + principalStretchDevPow[1][m] + principalStretchDevPow[2][m]) / 9.0);
                        }
                    }
                } else {
                    for (int m = 0; m < Nmax; m++) {
                        if (myMu[m] != 0) {
                            gab += myMu[m] * (-principalStretchDevPow[i][m] / 3.0 - principalStretchDevPow[n][m] / 3.0 + (principalStretchDevPow[0][m] + principalStretchDevPow[1][m] + principalStretchDevPow[2][m]) / 9.0);
                        }
                    }
                }
                // 2nd term in Eq. 2.68 of Simo and Taylor (1991)
                TensorUtils.addSymmetricTensorProduct(c, gab / J / 2.0, ma, mb);
            }
        }
    } else {
        // two are the same
        int i;
        if ((Math.abs(principalStretch2.get(0) - principalStretch2.get(1)) <= smallNum)) {
            i = 2;
        } else if (Math.abs(principalStretch2.get(0) - principalStretch2.get(2)) <= smallNum) {
            i = 1;
        } else {
            i = 0;
        }
        int j = (i + 1) % 3;
        int k = (j + 1) % 3;
        // D prime i - Eq. 2.46a of Simo and Taylor (1991)
        double Dpi = 8.0 * Math.pow(principalStretch.get(i), 3) - 2.0 * I1 * principalStretch.get(i) - 2.0 * I3 * Math.pow(principalStretch.get(i), -3);
        // D i - Eq. 2.16a of Simo and Taylor (1991)
        double Di = (principalStretch2.get(i) - principalStretch2.get(j)) * (principalStretch2.get(i) - principalStretch2.get(k));
        double beta3 = 0.0;
        double beta1 = 0.0;
        double g33 = 0.0;
        double g11 = 0.0;
        double g13 = 0.0;
        for (int n = 0; n < Nmax; n++) {
            if (myMu[n] != 0) {
                // beta_i in Eq. 2.65 of Simo and Taylor (1991)
                beta3 += myMu[n] / myAlpha[n] * (principalStretchDevPow[i][n] - (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 3.0);
                beta1 += myMu[n] / myAlpha[n] * (principalStretchDevPow[j][n] - (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 3.0);
                // gamma_AB in Eq. 2.66 of Simo and Taylor (1991)
                g33 += myMu[n] * (principalStretchDevPow[i][n] / 3.0 + (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 9.0);
                g11 += myMu[n] * (principalStretchDevPow[j][n] / 3.0 + (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 9.0);
                g13 += myMu[n] * (-principalStretchDevPow[i][n] / 3.0 - principalStretchDevPow[j][n] / 3.0 + (principalStretchDevPow[0][n] + principalStretchDevPow[1][n] + principalStretchDevPow[2][n]) / 9.0);
            }
        }
        // the matrix mi - Eq. 2.15 of Simo and Taylor (1991)
        ma.set(myB2);
        ma.scaledAdd(-principalStretch2.get(k), myB);
        ma.scaledAdd(-principalStretch2.get(j), myB);
        ma.scaledAdd(principalStretch2.get(j) * principalStretch2.get(k), SymmetricMatrix3d.IDENTITY);
        ma.scale(1.0 / Di);
        // Calculate dgm term in Eq. 2.48b and Eq 2.70 of Simo and Taylor (1991)
        TensorUtils.addTensorProduct4(c, 2.0 * (beta3 - beta1) / J / Di, myB);
        TensorUtils.addTensorProduct(c, -2.0 * (beta3 - beta1) / J / Di, myB);
        TensorUtils.addScaledIdentityProduct(c, 2.0 * (beta3 - beta1) / J * I3 / Di / principalStretch2.get(i));
        TensorUtils.addScaledIdentity(c, -2.0 * (beta3 - beta1) / J * I3 / Di / principalStretch2.get(i));
        TensorUtils.addSymmetricTensorProduct(c, 2.0 * (beta3 - beta1) / J / Di * principalStretch2.get(i), myB, ma);
        TensorUtils.addTensorProduct(c, -1.0 * (beta3 - beta1) / J / Di * Dpi * principalStretch.get(i), ma);
        TensorUtils.addSymmetricTensorProduct(c, -2.0 * (beta3 - beta1) / J / Di * I3 / principalStretch2.get(i), SymmetricMatrix3d.IDENTITY, ma);
        // Calculate other terms in Eq 2.70 of Simo and Taylor (1991)
        TensorUtils.addScaledIdentity(c, -2.0 * beta1 / J);
        myTmp.set(SymmetricMatrix3d.IDENTITY);
        myTmp.scaledAdd(-1.0, ma);
        TensorUtils.addTensorProduct(c, g11 / J, myTmp);
        TensorUtils.addTensorProduct(c, g33 / J, ma);
        TensorUtils.addSymmetricTensorProduct(c, g13 / J, ma, myTmp);
    }
    c.m00 += -p;
    c.m11 += -p;
    c.m22 += -p;
    c.m01 += p;
    c.m02 += p;
    c.m12 += p;
    c.m33 += -p;
    c.m44 += -p;
    c.m55 += -p;
    c.setLowerToUpper();
}

Example 17 with Matrix3d

use of maspack.matrix.Matrix3d in project artisynth_core by artisynth.

the class FungMaterial method main.

public static void main(String[] args) {
    FungMaterial mat = new FungMaterial();
    RotationMatrix3d R = new RotationMatrix3d();
    R.setRpy(1, 2, 3);
    SolidDeformation def = new SolidDeformation();
    Matrix3d Q = new Matrix3d();
    def.setF(new Matrix3d(1, 3, 5, 2, 1, 4, 6, 1, 2));
    Q.set(R);
    Matrix6d D = new Matrix6d();
    SymmetricMatrix3d sig = new SymmetricMatrix3d();
    // mat.computeStress (sig, pt, dt, null);
    // System.out.println ("sig=\n" + sig.toString ("%12.6f"));
    // mat.computeTangent (D, sig, pt, dt, null);
    // System.out.println ("D=\n" + D.toString ("%12.6f"));
    mat.computeStress(sig, def, Q, null);
    System.out.println("sig=\n" + sig.toString("%12.6f"));
    mat.computeTangent(D, sig, def, Q, null);
    System.out.println("D=\n" + D.toString("%12.6f"));
}

Example 18 with Matrix3d

use of maspack.matrix.Matrix3d in project artisynth_core by artisynth.

the class GenericMuscle method clone.

public GenericMuscle clone() {
    GenericMuscle mat = (GenericMuscle) super.clone();
    mat.myTmp = new Vector3d();
    mat.myMat = new Matrix3d();
    return mat;
}

Example 19 with Matrix3d

use of maspack.matrix.Matrix3d in project artisynth_core by artisynth.

the class MooneyRivlinMaterial method main.

public static void main(String[] args) {
    MooneyRivlinMaterial mat = new MooneyRivlinMaterial();
    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.setC10(1.2);
    mat.setC01(2.4);
    mat.setC11(3.5);
    mat.setC02(2.6);
    mat.setC20(1.9);
    mat.setBulkModulus(0.0);
    mat.computeStress(sig, def, Q, null);
    // pt.setStress (sig);
    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"));
}

Example 20 with Matrix3d

use of maspack.matrix.Matrix3d in project artisynth_core by artisynth.

the class BlemkerMuscle method clone.

public BlemkerMuscle clone() {
    BlemkerMuscle mat = (BlemkerMuscle) super.clone();
    mat.myTmp = new Vector3d();
    mat.myMat = new Matrix3d();
    return mat;
}