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

Example 16 with Matrix6d

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

the class CubicHyperelastic method main.

public static void main(String[] args) {
    CubicHyperelastic mat = new CubicHyperelastic();
    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.setG10(1.2);
    mat.setG30(3.5);
    mat.setG20(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 17 with Matrix6d

use of maspack.matrix.Matrix6d 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();
}

Example 18 with Matrix6d

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

the class FungMaterial method ddots.

public static Matrix6d ddots(Matrix6d a, Matrix6d b) {
    Matrix6d c = new Matrix6d();
    c.m00 = 2 * (a.m00 * b.m00 + a.m01 * b.m01 + a.m02 * b.m02 + 2 * a.m03 * b.m03 + 2 * a.m04 * b.m04 + 2 * a.m05 * b.m05);
    c.m01 = a.m00 * b.m01 + a.m11 * b.m01 + a.m01 * (b.m00 + b.m11) + a.m12 * b.m02 + a.m02 * b.m12 + 2 * a.m13 * b.m03 + 2 * a.m03 * b.m13 + 2 * a.m14 * b.m04 + 2 * a.m04 * b.m14 + 2 * a.m15 * b.m05 + 2 * a.m05 * b.m15;
    c.m02 = a.m12 * b.m01 + a.m00 * b.m02 + a.m22 * b.m02 + a.m01 * b.m12 + a.m02 * (b.m00 + b.m22) + 2 * a.m23 * b.m03 + 2 * a.m03 * b.m23 + 2 * a.m24 * b.m04 + 2 * a.m04 * b.m24 + 2 * a.m25 * b.m05 + 2 * a.m05 * b.m25;
    c.m03 = a.m13 * b.m01 + a.m23 * b.m02 + a.m00 * b.m03 + 2 * a.m33 * b.m03 + a.m01 * b.m13 + a.m02 * b.m23 + a.m03 * (b.m00 + 2 * b.m33) + 2 * a.m34 * b.m04 + 2 * a.m04 * b.m34 + 2 * a.m35 * b.m05 + 2 * a.m05 * b.m35;
    c.m04 = a.m14 * b.m01 + a.m24 * b.m02 + 2 * a.m34 * b.m03 + a.m00 * b.m04 + 2 * a.m44 * b.m04 + a.m01 * b.m14 + a.m02 * b.m24 + 2 * a.m03 * b.m34 + a.m04 * (b.m00 + 2 * b.m44) + 2 * a.m45 * b.m05 + 2 * a.m05 * b.m45;
    c.m05 = a.m15 * b.m01 + a.m25 * b.m02 + 2 * a.m35 * b.m03 + 2 * a.m45 * b.m04 + a.m00 * b.m05 + 2 * a.m55 * b.m05 + a.m01 * b.m15 + a.m02 * b.m25 + 2 * a.m03 * b.m35 + 2 * a.m04 * b.m45 + a.m05 * (b.m00 + 2 * b.m55);
    c.m11 = 2 * (a.m01 * b.m01 + a.m11 * b.m11 + a.m12 * b.m12 + 2 * a.m13 * b.m13 + 2 * a.m14 * b.m14 + 2 * a.m15 * b.m15);
    c.m12 = a.m02 * b.m01 + a.m01 * b.m02 + a.m11 * b.m12 + a.m22 * b.m12 + a.m12 * (b.m11 + b.m22) + 2 * a.m23 * b.m13 + 2 * a.m13 * b.m23 + 2 * a.m24 * b.m14 + 2 * a.m14 * b.m24 + 2 * a.m25 * b.m15 + 2 * a.m15 * b.m25;
    c.m13 = a.m03 * b.m01 + a.m23 * b.m12 + a.m01 * b.m03 + a.m11 * b.m13 + 2 * a.m33 * b.m13 + a.m12 * b.m23 + a.m13 * (b.m11 + 2 * b.m33) + 2 * a.m34 * b.m14 + 2 * a.m14 * b.m34 + 2 * a.m35 * b.m15 + 2 * a.m15 * b.m35;
    c.m14 = a.m04 * b.m01 + a.m24 * b.m12 + 2 * a.m34 * b.m13 + a.m01 * b.m04 + a.m11 * b.m14 + 2 * a.m44 * b.m14 + a.m12 * b.m24 + 2 * a.m13 * b.m34 + a.m14 * (b.m11 + 2 * b.m44) + 2 * a.m45 * b.m15 + 2 * a.m15 * b.m45;
    c.m15 = a.m05 * b.m01 + a.m25 * b.m12 + 2 * a.m35 * b.m13 + 2 * a.m45 * b.m14 + a.m01 * b.m05 + a.m11 * b.m15 + 2 * a.m55 * b.m15 + a.m12 * b.m25 + 2 * a.m13 * b.m35 + 2 * a.m14 * b.m45 + a.m15 * (b.m11 + 2 * b.m55);
    c.m22 = 2 * (a.m02 * b.m02 + a.m12 * b.m12 + a.m22 * b.m22 + 2 * a.m23 * b.m23 + 2 * a.m24 * b.m24 + 2 * a.m25 * b.m25);
    c.m23 = a.m03 * b.m02 + a.m13 * b.m12 + a.m23 * b.m22 + a.m02 * b.m03 + a.m12 * b.m13 + a.m22 * b.m23 + 2 * a.m33 * b.m23 + 2 * a.m23 * b.m33 + 2 * a.m34 * b.m24 + 2 * a.m24 * b.m34 + 2 * a.m35 * b.m25 + 2 * a.m25 * b.m35;
    c.m24 = a.m04 * b.m02 + a.m14 * b.m12 + a.m24 * b.m22 + 2 * a.m34 * b.m23 + a.m02 * b.m04 + a.m12 * b.m14 + a.m22 * b.m24 + 2 * a.m44 * b.m24 + 2 * a.m23 * b.m34 + 2 * a.m24 * b.m44 + 2 * a.m45 * b.m25 + 2 * a.m25 * b.m45;
    c.m25 = a.m05 * b.m02 + a.m15 * b.m12 + a.m25 * b.m22 + 2 * a.m35 * b.m23 + 2 * a.m45 * b.m24 + a.m02 * b.m05 + a.m12 * b.m15 + a.m22 * b.m25 + 2 * a.m55 * b.m25 + 2 * a.m23 * b.m35 + 2 * a.m24 * b.m45 + 2 * a.m25 * b.m55;
    c.m33 = 2 * (a.m03 * b.m03 + a.m13 * b.m13 + a.m23 * b.m23 + 2 * a.m33 * b.m33 + 2 * a.m34 * b.m34 + 2 * a.m35 * b.m35);
    c.m34 = a.m04 * b.m03 + a.m14 * b.m13 + a.m24 * b.m23 + 2 * a.m34 * b.m33 + a.m03 * b.m04 + a.m13 * b.m14 + a.m23 * b.m24 + 2 * a.m33 * b.m34 + 2 * a.m44 * b.m34 + 2 * a.m34 * b.m44 + 2 * a.m45 * b.m35 + 2 * a.m35 * b.m45;
    c.m35 = a.m05 * b.m03 + a.m15 * b.m13 + a.m25 * b.m23 + 2 * a.m35 * b.m33 + 2 * a.m45 * b.m34 + a.m03 * b.m05 + a.m13 * b.m15 + a.m23 * b.m25 + 2 * a.m33 * b.m35 + 2 * a.m55 * b.m35 + 2 * a.m34 * b.m45 + 2 * a.m35 * b.m55;
    c.m44 = 2 * (a.m04 * b.m04 + a.m14 * b.m14 + a.m24 * b.m24 + 2 * a.m34 * b.m34 + 2 * a.m44 * b.m44 + 2 * a.m45 * b.m45);
    c.m45 = a.m05 * b.m04 + a.m15 * b.m14 + a.m25 * b.m24 + 2 * a.m35 * b.m34 + 2 * a.m45 * b.m44 + a.m04 * b.m05 + a.m14 * b.m15 + a.m24 * b.m25 + 2 * a.m34 * b.m35 + 2 * a.m44 * b.m45 + 2 * a.m55 * b.m45 + 2 * a.m45 * b.m55;
    c.m55 = 2 * (a.m05 * b.m05 + a.m15 * b.m15 + a.m25 * b.m25 + 2 * a.m35 * b.m35 + 2 * a.m45 * b.m45 + 2 * a.m55 * b.m55);
    return c;
}

Example 19 with Matrix6d

use of maspack.matrix.Matrix6d 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"));
}

Example 20 with Matrix6d

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

the class AnisotropicLinearMaterial method createIsotropicStiffness.

public static Matrix6d createIsotropicStiffness(double E, double nu) {
    Matrix6d D = new Matrix6d();
    double lam = E * nu / ((1 + nu) * (1 - 2 * nu));
    double mu = E / (2 * (1 + nu));
    D.m00 = lam + 2 * mu;
    D.m01 = lam;
    D.m02 = lam;
    D.m10 = lam;
    D.m11 = lam + 2 * mu;
    D.m12 = lam;
    D.m20 = lam;
    D.m21 = lam;
    D.m22 = lam + 2 * mu;
    D.m33 = mu;
    D.m44 = mu;
    D.m55 = mu;
    return D;
}