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

Example 26 with SymmetricMatrix3d

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

the class PointDistributor method getPrincipalAxes.

public static RigidTransform3d getPrincipalAxes(PolygonalMesh mesh) {
    Vector3d mov1 = new Vector3d();
    Vector3d mov2 = new Vector3d();
    Vector3d pov = new Vector3d();
    double vol = mesh.computeVolumeIntegrals(mov1, mov2, pov);
    double mass = vol;
    Point3d cov = new Point3d();
    // center of volume
    cov.scale(1.0 / vol, mov1);
    // [c], skew symmetric
    Matrix3d covMatrix = new Matrix3d(0, -cov.z, cov.y, cov.z, 0, -cov.x, -cov.y, cov.x, 0);
    // J
    Matrix3d J = new Matrix3d((mov2.y + mov2.z), -pov.z, -pov.y, -pov.z, (mov2.x + mov2.z), -pov.x, -pov.y, -pov.x, (mov2.x + mov2.y));
    // Jc = J + m[c][c]
    Matrix3d Jc = new Matrix3d();
    Jc.mul(covMatrix, covMatrix);
    Jc.scale(mass);
    Jc.add(J);
    // Compute eigenvectors and eigenvlaues of Jc
    SymmetricMatrix3d JcSymmetric = new SymmetricMatrix3d(Jc);
    Vector3d lambda = new Vector3d();
    Matrix3d U = new Matrix3d();
    JcSymmetric.getEigenValues(lambda, U);
    // Construct the rotation matrix
    RotationMatrix3d R = new RotationMatrix3d();
    R.set(U);
    lambda.absolute();
    if (lambda.x > lambda.y && lambda.z > lambda.y) {
        R.rotateZDirection(new Vector3d(R.m01, R.m11, R.m21));
    } else if (lambda.x > lambda.z && lambda.y > lambda.z) {
        R.rotateZDirection(new Vector3d(R.m00, R.m10, R.m20));
    }
    return (new RigidTransform3d(cov, R));
}

Example 27 with SymmetricMatrix3d

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

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

the class FullBlemkerMuscle method computeTangent.

public void computeTangent(Matrix6d D, SymmetricMatrix3d stress, double excitation, Vector3d dir0, SolidDeformation def, FemMaterial baseMat) {
    Vector3d a = myTmp;
    def.getF().mul(a, dir0);
    double mag = a.norm();
    a.scale(1 / mag);
    double J = def.getDetF();
    double Ji = 1 / J;
    double lamd = mag * Math.pow(J, -1.0 / 3.0);
    double I4 = lamd * lamd;
    // BEGIN G1, G2
    double I1 = 0;
    double I2 = 0;
    double I5 = 0;
    // calculate deviatoric left Cauchy-Green tensor
    def.computeDevLeftCauchyGreen(myB);
    // calculate square of B
    myB2.mulTransposeLeft(myB);
    Vector3d Ba = new Vector3d();
    myB.mul(Ba, a);
    // Invariants of deviatoric part of B
    I1 = myB.trace();
    I2 = 0.5 * (I1 * I1 - myB2.trace());
    I5 = I4 * Ba.dot(a);
    // calculate new invariants
    double g = I5 / (I4 * I4) - 1;
    double b1 = (g > 0 ? Math.sqrt(g) : 0);
    double b2 = acosh(0.5 * (I1 * I4 - I5) / lamd);
    // calculate omage (w)
    double w = 0.5 * (I1 * I4 - I5) / lamd;
    // set beta and ksi to their limit values
    double beta = 1.0;
    double ksi = -1.0 / 3.0;
    // if w not equals unity, we can calculate beta and ksi
    if (w > 1.0001) {
        beta = b2 / Math.sqrt(w * w - 1);
        ksi = (1.0 / (w * w - 1)) * (1 - w * b2 / Math.sqrt(w * w - 1));
    }
    // -- A. matrix contribution --
    // calculate derivatives for F1
    double F1D4 = -2 * myG1 * (I5 / (I4 * I4 * I4));
    double F1D5 = myG1 / (I4 * I4);
    double F1D44 = 6 * myG1 * (I5 / (I4 * I4 * I4 * I4));
    double F1D45 = -2 * myG1 / (I4 * I4 * I4);
    // calculate derivatives for F2
    double F2D1 = myG2 * beta * lamd;
    double F2D4 = myG2 * beta * (I1 * I4 + I5) * 0.5 * Math.pow(I4, -1.5);
    double F2D5 = -myG2 * beta / lamd;
    double F2D11 = ksi * myG2 * I4 * 0.5;
    double F2D44 = 2.0 * myG2 * ksi * Math.pow(0.25 * (I1 * I4 + I5) / Math.pow(I4, 1.5), 2) - myG2 * beta * (0.25 * (I1 * I4 + 3 * I5) / Math.pow(I4, 2.5));
    double F2D55 = 0.5 * myG2 * ksi / I4;
    double F2D14 = myG2 * beta * 0.5 / lamd + myG2 * ksi * (I1 * I4 + I5) * 0.25 / I4;
    double F2D15 = -0.5 * myG2 * ksi;
    double F2D45 = myG2 * beta * 0.5 * Math.pow(I4, -1.5) - myG2 * ksi * 0.25 * (I1 * I4 + I5) / (I4 * I4);
    // calculate derivatives for F3
    // these terms are proposed to fix the zero-stress problem
    // double F3D4  = 9.0*myG3*0.125*log(I4)/I4;
    // double F3D44 = 9.0*myG3*0.125*(1 - log(I4))/(I4*I4);
    // END G1, G2
    double P1 = myExpStressCoeff;
    double P2 = myUncrimpingFactor;
    double lofl = myOptLambda;
    double Fa = 0, Fp;
    double FaDl = 0, FpDl;
    if (!myP3P4Valid) {
        myP3 = P1 * P2 * Math.exp(P2 * (myMaxLambda / lofl - 1));
        myP4 = P1 * (Math.exp(P2 * (myMaxLambda / lofl - 1)) - 1) - myP3 * myMaxLambda / lofl;
        myP3P4Valid = true;
    }
    // calculate passive fiber force
    if (myZeroForceBelowLamOptP && lamd <= lofl) {
        // Fp is zero if less than lofl
        Fp = 0;
        FpDl = 0;
    } else if (lamd < myMaxLambda) {
        Fp = P1 * (Math.exp(P2 * (lamd / lofl - 1)) - 1);
        FpDl = P1 * P2 * Math.exp(P2 * (lamd / lofl - 1)) / lofl;
    } else {
        Fp = myP3 * lamd / lofl + myP4;
        FpDl = myP3 / lofl;
    }
    // calculate active fiber force
    if ((lamd <= 0.4 * lofl) || (lamd >= 1.6 * lofl)) {
        // FEBio has added this part to make sure that
        // Fa is zero outside the range [0.4, 1.6] *lofl
        Fa = 0;
        FaDl = 0;
    } else {
        if (lamd <= 0.6 * lofl) {
            Fa = 9 * square(lamd / lofl - 0.4);
            FaDl = 18 * (lamd / lofl - 0.4) / lofl;
        } else if (lamd >= 1.4 * lofl) {
            Fa = 9 * square(lamd / lofl - 1.6);
            FaDl = 18 * (lamd / lofl - 1.6) / lofl;
        } else if ((lamd >= 0.6 * lofl) && (lamd <= 1.4 * lofl)) {
            Fa = 1 - 4 * square(1 - lamd / lofl);
            FaDl = 8 * (1 - lamd / lofl) / lofl;
        }
    }
    // calculate total fiber force
    double FfDl = myMaxStress * (Fp + excitation * Fa) / lofl;
    double FfD4 = 0.5 * FfDl / lamd;
    double FfDll = myMaxStress * (FpDl + excitation * FaDl) / lofl;
    double FfD44 = 0.25 * (FfDll - FfDl / lamd) / I4;
    double W1, W2, W4, W5;
    // add all derivatives
    W1 = F2D1;
    W2 = 0;
    // + F3D4
    W4 = F1D4 + F2D4 + FfD4;
    W5 = F1D5 + F2D5;
    // calculate second derivatives
    double W11, W12, W22, W14, W24, W15, W25, W44, W45, W55;
    W11 = F2D11;
    W12 = 0;
    W22 = 0;
    W14 = F2D14;
    W24 = 0;
    W15 = F2D15;
    W25 = 0;
    // + F3D44
    W44 = F1D44 + F2D44 + FfD44;
    W45 = F1D45 + F2D45;
    W55 = F2D55;
    // System.out.printf ("W1=%g W2=%g W4=%g W5=%g\n", W1, W2, W4, W5);
    // System.out.printf ("W11=%g W12=%g W22=%g W14=%g W24=%g\n",
    // W11, W12, W22, W14, W24);
    // System.out.printf ("W15=%g W25=%g W44=%g W45=%g W55=%g\n",
    // W15, W25, W44, W45, W55);
    // calculate dWdC:C
    double WCC = W1 * I1 + 2 * W2 * I2 + W4 * I4 + 2 * W5 * I5;
    // calculate C:d2WdCdC:C
    double CW2CCC = (W11 * I1 + W12 * I1 * I1 + W2 * I1 + 2 * W12 * I2 + 2 * W22 * I1 * I2 + W14 * I4 + W24 * I1 * I4 + 2 * W15 * I5 + 2 * W25 * I1 * I5) * I1 - (W12 * I1 + 2 * W22 * I2 + W2 + W24 * I4 + 2 * W25 * I5) * (I1 * I1 - 2 * I2) + (W14 * I1 + 2 * W24 * I2 + W44 * I4 + 2 * W45 * I5) * I4 + (W15 * I1 + 2 * W25 * I2 + W45 * I4 + 2 * W55 * I5) * 2 * I5 + 2 * W5 * I5;
    SymmetricMatrix3d AA = new SymmetricMatrix3d();
    SymmetricMatrix3d AB = new SymmetricMatrix3d();
    SymmetricMatrix3d WCCC = myMat;
    AA.dyad(a);
    AB.symmetricDyad(a, Ba);
    WCCC.scale(W11 * I1 + W12 * I1 * I1 + W2 * I1 + 2 * W12 * I2 + 2 * W22 * I1 * I2 + W14 * I4 + W24 * I1 * I4 + 2 * W15 * I5 + 2 * W25 * I1 * I5, myB);
    WCCC.scaledAdd(-(W12 * I1 + 2 * W22 * I2 + W2 + W24 * I4 + 2 * W25 * I5), myB2, WCCC);
    WCCC.scaledAdd((W14 * I1 + 2 * W24 * I2 + W44 * I4 + 2 * W45 * I5) * I4, AA, WCCC);
    WCCC.scaledAdd((W15 * I1 + 2 * W25 * I2 + W45 * I4 + 2 * W55 * I5 + W5) * I4, AB, WCCC);
    D.setZero();
    TensorUtils.addTensorProduct(D, (W11 + 2.0 * W12 * I1 + W2 + W22 * I1 * I1) * 4 * Ji, myB);
    TensorUtils.addSymmetricTensorProduct(D, -(W12 + W22 * I1) * 4 * Ji, myB, myB2);
    TensorUtils.addTensorProduct(D, W22 * 4 * Ji, myB2);
    TensorUtils.addSymmetricTensorProduct(D, (W14 + W24 * I1) * I4 * 4 * Ji, myB, AA);
    TensorUtils.addSymmetricTensorProduct(D, (W15 + W25 * I1) * I4 * 4 * Ji, myB, AB);
    TensorUtils.addSymmetricTensorProduct(D, (-W24 * I4) * 4 * Ji, myB2, AA);
    TensorUtils.addTensorProduct(D, (W44 * I4 * I4) * 4 * Ji, AA);
    TensorUtils.addSymmetricTensorProduct(D, (W45 * I4) * I4 * 4 * Ji, AA, AB);
    TensorUtils.addTensorProduct(D, (W55) * I4 * I4 * 4 * Ji, AB);
    TensorUtils.addSymmetricTensorProduct4(D, W5 * I4 * 4 * Ji, AA, myB);
    TensorUtils.addScaledIdentityProduct(D, 4 / 9.0 * Ji * (CW2CCC - WCC));
    WCCC.scale(-4 / 3.0 * Ji);
    TensorUtils.addSymmetricIdentityProduct(D, WCCC);
    TensorUtils.addScaledIdentity(D, 4 / 3.0 * Ji * WCC);
    // compute stress (in myMat) due to this material
    myMat.scale(W1 + W2 * I1, myB);
    myMat.scaledAdd(-W2, myB2, myMat);
    myMat.scaledAdd(I4 * W4, AA);
    myMat.scaledAdd(I4 * W5, AB);
    myMat.deviator();
    myMat.scale(2.0 / J);
    // myMat.set (def.getStrain());
    // myMat.deviator();
    myMat.scale(-2.0 / 3.0);
    TensorUtils.addSymmetricIdentityProduct(D, myMat);
    D.setLowerToUpper();
}

Example 29 with SymmetricMatrix3d

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

the class FungMaterial method computeStress.

public void computeStress(SymmetricMatrix3d sigma, SolidDeformation def, Matrix3d Q, FemMaterial baseMat) {
    sigma.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();
    SymmetricMatrix3d sigmaFung = new SymmetricMatrix3d();
    // 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 stress
    SymmetricMatrix3d bmi = new SymmetricMatrix3d();
    bmi.sub(myB, SymmetricMatrix3d.IDENTITY);
    for (int i = 0; i < 3; i++) {
        // s += 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);
        sigmaFung.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.;
            sigmaFung.scaledAdd(lam[i][j] / 2.0 * (K[i] - 1.0) * K[j], A[j]);
            sigmaFung.scaledAdd(lam[i][j] / 2.0 * (K[j] - 1.0) * K[i], A[i]);
        }
    }
    sigmaFung.scale(eQ / (2.0 * J));
    sigmaFung.deviator();
    sigmaFung.m10 = sigmaFung.m01;
    sigmaFung.m20 = sigmaFung.m02;
    sigmaFung.m21 = sigmaFung.m12;
    sigma.add(sigmaFung);
    sigma.m00 += avgp;
    sigma.m11 += avgp;
    sigma.m22 += avgp;
}

Example 30 with SymmetricMatrix3d

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

the class FungMaterial method clone.

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