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

Example 6 with SymmetricMatrix3d

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

the class OgdenMaterial method main.

public static void main(String[] args) {
    OgdenMaterial mat = new OgdenMaterial();
    SolidDeformation def = new SolidDeformation();
    Matrix3d Q = new Matrix3d();
    def.setF(new Matrix3d(1.1, 0.1, 0.2, 0.3, 0.8, 0.23, 0.3, 0.1, 1.5));
    Matrix6d D = new Matrix6d();
    SymmetricMatrix3d sig = new SymmetricMatrix3d();
    // double[] alpha = {2.0, 2.0, 2.0, 2.0, 2.0, 2.0};
    // double[] mu = {200.0, 0.0, 0.0, 0.0, 0.0, 0.0};
    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 7 with SymmetricMatrix3d

use of maspack.matrix.SymmetricMatrix3d 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 8 with SymmetricMatrix3d

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

the class QLVBehavior method computeStress.

public void computeStress(SymmetricMatrix3d sigma, ViscoelasticState state) {
    QLVState qlvstate = (QLVState) state;
    double h = qlvstate.myH;
    if (h == 0) {
        return;
    }
    double[] GHPrev = qlvstate.myGHPrev;
    double[] S = qlvstate.myS;
    SymmetricMatrix3d deltaSigma = new SymmetricMatrix3d();
    deltaSigma.sub(sigma, qlvstate.mySigmaPrev);
    qlvstate.mySigmaSave.set(sigma);
    sigma.scale(myGamma0);
    for (int i = 0; i < N_MAX; i++) {
        double H00 = S[i] * deltaSigma.m00 + GHPrev[i * 6];
        double H11 = S[i] * deltaSigma.m11 + GHPrev[i * 6 + 1];
        double H22 = S[i] * deltaSigma.m22 + GHPrev[i * 6 + 2];
        double H01 = S[i] * deltaSigma.m01 + GHPrev[i * 6 + 3];
        double H02 = S[i] * deltaSigma.m02 + GHPrev[i * 6 + 4];
        double H12 = S[i] * deltaSigma.m12 + GHPrev[i * 6 + 5];
        double gamma = myGamma[i];
        sigma.m00 += gamma * H00;
        sigma.m11 += gamma * H11;
        sigma.m22 += gamma * H22;
        sigma.m01 += gamma * H01;
        sigma.m02 += gamma * H02;
        sigma.m12 += gamma * H12;
    }
    sigma.m10 = sigma.m01;
    sigma.m20 = sigma.m02;
    sigma.m21 = sigma.m12;
}

Example 9 with SymmetricMatrix3d

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

the class QLVBehavior method advanceState.

public void advanceState(ViscoelasticState state, double t0, double t1) {
    QLVState qlvstate = (QLVState) state;
    double h = t1 - t0;
    qlvstate.myH = h;
    SymmetricMatrix3d deltaSigma = new SymmetricMatrix3d();
    deltaSigma.sub(qlvstate.mySigmaSave, qlvstate.mySigmaPrev);
    if (t0 >= 0) {
        qlvstate.mySigmaPrev.set(qlvstate.mySigmaSave);
    }
    double[] S = qlvstate.myS;
    double[] GHPrev = qlvstate.myGHPrev;
    myTangentScale = myGamma0;
    for (int i = 0; i < N_MAX; i++) {
        double g = Math.exp(-h / myTau[i]);
        double gH;
        int idx = i * 6;
        gH = g * (S[i] * deltaSigma.m00 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        gH = g * (S[i] * deltaSigma.m11 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        gH = g * (S[i] * deltaSigma.m22 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        gH = g * (S[i] * deltaSigma.m01 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        gH = g * (S[i] * deltaSigma.m02 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        gH = g * (S[i] * deltaSigma.m12 + GHPrev[idx]);
        GHPrev[idx++] = gH;
        S[i] = (1.0 - g) / (h / myTau[i]);
        myTangentScale += myGamma[i] * S[i];
    }
}

Example 10 with SymmetricMatrix3d

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

the class CubicHyperelastic method clone.

public CubicHyperelastic clone() {
    CubicHyperelastic mat = (CubicHyperelastic) super.clone();
    mat.myB = new SymmetricMatrix3d();
    mat.myTmp = new SymmetricMatrix3d();
    return mat;
}