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

Example 56 with Matrix3d

use of maspack.matrix.Matrix3d 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 57 with Matrix3d

use of maspack.matrix.Matrix3d 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 58 with Matrix3d

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

Example 59 with Matrix3d

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

the class OgdenMaterial method computeStress.

public void computeStress(SymmetricMatrix3d sigma, SolidDeformation def, Matrix3d Q, FemMaterial baseMat) {
    double J = def.getDetF();
    double avgp = def.getAveragePressure();
    sigma.setZero();
    // Calculate Deviatoric left Cauchy-Green tensor
    def.computeDevLeftCauchyGreen(myB);
    Vector3d principalStretch = new Vector3d();
    Vector3d principalStretch2 = 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 principal stresses
    for (int i = 0; i < 3; i++) {
        for (int n = 0; n < Nmax; n++) {
            if (myMu[n] != 0) {
                sigma.set(i, i, sigma.get(i, i) + myMu[n] / myAlpha[n] / J * (Math.pow(principalStretch.get(i), myAlpha[n])));
            }
        }
    }
    // Calculate stress tensor from principal stresses and directions
    sigma.mulLeftAndTransposeRight(principalDirection);
    sigma.deviator();
    sigma.m00 += avgp;
    sigma.m11 += avgp;
    sigma.m22 += avgp;
}

Example 60 with Matrix3d

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

the class RotAxisFrameMaterial method computeDFdq.

public void computeDFdq(Matrix6d Jq, RigidTransform3d X21, Twist vel21, RigidTransform3d initialX21, boolean symmetric) {
    AxisAngle axisAng = new AxisAngle();
    X21.R.getAxisAngle(axisAng);
    Matrix3d U = new Matrix3d();
    Jq.setZero();
    computeU(U, axisAng, symmetric);
    Jq.m00 = myK.x;
    Jq.m11 = myK.y;
    Jq.m22 = myK.z;
    // set lower right submatrix to myRotK*U
    Jq.m33 = myRotK * U.m00;
    Jq.m34 = myRotK * U.m01;
    Jq.m35 = myRotK * U.m02;
    Jq.m43 = myRotK * U.m10;
    Jq.m44 = myRotK * U.m11;
    Jq.m45 = myRotK * U.m12;
    Jq.m53 = myRotK * U.m20;
    Jq.m54 = myRotK * U.m21;
    Jq.m55 = myRotK * U.m22;
}