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

Example 21 with SymmetricMatrix3d

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

the class AuxMaterialElementDesc method computeStress.

@Override
public void computeStress(SymmetricMatrix3d sigma, SolidDeformation def, IntegrationPoint3d pt, IntegrationData3d dt, FemMaterial baseMat) {
    FemMaterial mat = getEffectiveMaterial();
    if (mat != null) {
        double frac = myFrac;
        if (myFracs != null) {
            frac = myFracs[pt.getNumber()];
        }
        if (frac > 0) {
            Matrix3d Q = (dt.myFrame != null ? dt.myFrame : Matrix3d.IDENTITY);
            mat.computeStress(sigma, def, Q, baseMat);
            sigma.scale(frac);
        } else {
            sigma.setZero();
        }
    }
}

Example 22 with SymmetricMatrix3d

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

the class PolygonalMesh method computeInertia.

/**
 * Computes a spatial inertia for volume defined by this mesh, assuming a
 * constant density. It is assumed that the mesh is closed, although small
 * holes in the mesh should not affect the calculation that much.  Is is
 * also assumed that the faces are oriented counter-clockwise about their
 * outward-facing normal.
 *
 * The code for this was taken from vclip, by Brian Mirtich. See "Fast and
 * Accurate Computation of Polyhedral Mass Properties," Brian Mirtich,
 * journal of graphics tools, volume 1, number 2, 1996.
 *
 * @param M returns the computed spatial inertia
 * @param density
 * density of the volume
 * @return the volume of the mesh
 */
public double computeInertia(SpatialInertia M, double density) {
    Vector3d mov1 = new Vector3d();
    Vector3d mov2 = new Vector3d();
    Vector3d pov = new Vector3d();
    double vol = computeVolumeIntegrals(mov1, mov2, pov);
    Point3d cov = new Point3d();
    // center of volume
    cov.scale(1.0 / vol, mov1);
    double mass = density * vol;
    SymmetricMatrix3d J = new SymmetricMatrix3d(mov2.y + mov2.z, mov2.x + mov2.z, mov2.x + mov2.y, -pov.z, -pov.y, -pov.x);
    J.scale(density);
    // J contains the mass[com][com] term; remove this:
    J.m00 -= mass * (cov.z * cov.z + cov.y * cov.y);
    J.m11 -= mass * (cov.z * cov.z + cov.x * cov.x);
    J.m22 -= mass * (cov.y * cov.y + cov.x * cov.x);
    J.m01 += mass * cov.x * cov.y;
    J.m02 += mass * cov.x * cov.z;
    J.m12 += mass * cov.z * cov.y;
    J.m10 = J.m01;
    J.m20 = J.m02;
    J.m21 = J.m12;
    M.set(mass, J, cov);
    return vol;
}

Example 23 with SymmetricMatrix3d

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

the class FemModel3d method getMassMatrixValues.

@Override
public void getMassMatrixValues(SparseBlockMatrix M, VectorNd f, double t) {
    int bk;
    for (int k = 0; k < myNodes.size(); k++) {
        FemNode3d n = myNodes.get(k);
        if ((bk = n.getSolveIndex()) != -1) {
            n.getEffectiveMass(M.getBlock(bk, bk), t);
            n.getEffectiveMassForces(f, t, M.getBlockRowOffset(bk));
        }
    }
    if (myFrameRelativeP) {
        // angular velocity in body coords
        Vector3d wbod = new Vector3d();
        RotationMatrix3d R = myFrame.getPose().R;
        double[] fbuf = f.getBuffer();
        wbod.inverseTransform(R, myFrame.getVelocity().w);
        Vector3d w = myFrame.getVelocity().w;
        // update inertia in Frame
        double mass = 0;
        Point3d com = new Point3d();
        SymmetricMatrix3d J = new SymmetricMatrix3d();
        // extra fictitious forces for spatial inertia due to nodal velocity
        Vector3d fv = new Vector3d();
        Vector3d fw = new Vector3d();
        Vector3d fn = new Vector3d();
        Vector3d tmp = new Vector3d();
        // local node position rotated to world
        Point3d c = new Point3d();
        // local node velocity rotated to world
        Vector3d v = new Vector3d();
        int bf = myFrame.getSolveIndex();
        for (int k = 0; k < myNodes.size(); k++) {
            FemNode3d n = myNodes.get(k);
            if ((bk = n.getSolveIndex()) != -1) {
                c.transform(R, n.getLocalPosition());
                v.transform(R, n.getLocalVelocity());
                double m = n.getEffectiveMass();
                // System.out.println ("m " + k + " " + m);
                com.scaledAdd(m, c);
                mass += m;
                SpatialInertia.addPointRotationalInertia(J, m, c);
                if (useFrameRelativeCouplingMasses) {
                    Matrix3x6Block blk = (Matrix3x6Block) M.getBlock(bk, bf);
                    Matrix6x3Block blkT = (Matrix6x3Block) M.getBlock(bf, bk);
                    setCouplingMass(blk, m, R, c);
                    blkT.transpose(blk);
                    // compute fictitious forces for nodes, and accumulate nodal
                    // velocity fictitious force terms for spatial inertia
                    // tmp = 2*m (w X v)
                    tmp.cross(w, v);
                    tmp.scale(2 * m);
                    // fv += 2*m (w X v)
                    fv.add(tmp);
                    // fw += 2*m (c X w X v)
                    fw.crossAdd(c, tmp, fw);
                    // fn = 2*m (w X v)
                    fn.set(tmp);
                    // tmp = m*w X c
                    tmp.cross(w, c);
                    tmp.scale(m);
                    // fn += m (w X w X c)
                    fn.crossAdd(w, tmp, fn);
                    fn.inverseTransform(R);
                    // set fictitious force terms for node
                    int idx = M.getBlockRowOffset(bk);
                    fbuf[idx++] = -fn.x;
                    fbuf[idx++] = -fn.y;
                    fbuf[idx++] = -fn.z;
                }
            }
        }
        com.scale(1 / mass);
        SpatialInertia.addPointRotationalInertia(J, -mass, com);
        SpatialInertia S = new SpatialInertia();
        S.set(mass, J, com);
        // Frame keeps inertia in local coords
        S.inverseTransform(R);
        myFrame.setInertia(S);
        myFrame.getMass(M.getBlock(bf, bf), t);
        myFrame.getEffectiveMassForces(f, t, M.getBlockRowOffset(bf));
        if (useFrameRelativeCouplingMasses) {
            int idx = M.getBlockRowOffset(bf);
            fbuf[idx++] -= fv.x;
            fbuf[idx++] -= fv.y;
            fbuf[idx++] -= fv.z;
            fbuf[idx++] -= fw.x;
            fbuf[idx++] -= fw.y;
            fbuf[idx++] -= fw.z;
        }
        if (frameMassFraction != 0) {
            scaleFrameNodeMasses(M, f, 1.0 - frameMassFraction);
        }
    }
}

Example 24 with SymmetricMatrix3d

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

the class FemModel3d method computeStressAndStiffness.

// DIVBLK
private void computeStressAndStiffness(FemElement3d e, FemMaterial mat, Matrix6d D, IncompMethod softIncomp) {
    IntegrationPoint3d[] ipnts = e.getIntegrationPoints();
    IntegrationData3d[] idata = e.getIntegrationData();
    FemNode3d[] nodes = e.getNodes();
    if (D != null) {
        D.setZero();
    }
    SolidDeformation def = new SolidDeformation();
    // ===========================================
    // linear material optimizations
    // ===========================================
    // potentially update cached linear material
    // internally updates
    StiffnessWarper3d warper = e.getStiffnessWarper();
    // if there is cached linear material, then apply
    if (!warper.isCacheEmpty()) {
        // compute warping rotation
        warper.computeWarpingRotation(e);
        // add force and stiffness
        for (int i = 0; i < nodes.length; i++) {
            int bi = nodes[i].getSolveIndex();
            if (bi != -1) {
                FemNode3d n = nodes[i];
                if (!myStiffnessesValidP) {
                    for (int j = 0; j < nodes.length; j++) {
                        int bj = nodes[j].getSolveIndex();
                        if (!mySolveMatrixSymmetricP || bj >= bi) {
                            warper.addNodeStiffness(e.myNbrs[i][j].getK(), i, j);
                        }
                    }
                }
                // add node force
                warper.addNodeForce(n.myInternalForce, i, nodes);
            }
        }
        if (myComputeNodalStress || (myComputeNodalStrain && mat.isLinear())) {
            // estimate at warping point
            RotationMatrix3d R = warper.getRotation();
            IntegrationPoint3d wpnt = e.getWarpingPoint();
            IntegrationData3d wdata = e.getWarpingData();
            wpnt.computeJacobianAndGradient(nodes, wdata.myInvJ0);
            def.clear();
            def.setF(wpnt.F);
            def.setAveragePressure(0);
            def.setR(R);
            SymmetricMatrix3d sigma = new SymmetricMatrix3d();
            SymmetricMatrix3d tmp = new SymmetricMatrix3d();
            // compute nodal stress at wpnt
            if (myComputeNodalStress) {
                // compute linear stress
                mat.computeStress(tmp, def, null, null);
                sigma.add(tmp);
                for (AuxiliaryMaterial amat : e.getAuxiliaryMaterials()) {
                    amat.computeStress(tmp, def, wpnt, e.getWarpingData(), mat);
                    sigma.add(tmp);
                }
                // distribute stress to nodes
                for (int i = 0; i < nodes.length; i++) {
                    nodes[i].addScaledStress(1.0 / nodes[i].numAdjacentElements(), sigma);
                }
            }
            if (myComputeNodalStrain && mat.isLinear()) {
                // Cauchy strain at warping point
                if (mat.isCorotated()) {
                    // remove rotation from F
                    sigma.mulTransposeLeftSymmetric(R, wpnt.F);
                } else {
                    sigma.setSymmetric(wpnt.F);
                }
                sigma.m00 -= 1;
                sigma.m11 -= 1;
                sigma.m22 -= 1;
                // distribute strain to nodes
                for (int i = 0; i < nodes.length; i++) {
                    nodes[i].addScaledStrain(1.0 / nodes[i].numAdjacentElements(), sigma);
                }
            }
        }
    // stress or strain
    }
    // exit early if no non-linear materials
    boolean linearOnly = mat.isLinear();
    if (linearOnly) {
        for (AuxiliaryMaterial amat : e.getAuxiliaryMaterials()) {
            if (!amat.isLinear()) {
                linearOnly = false;
                break;
            }
        }
    }
    if (linearOnly) {
        return;
    }
    // we have some non-linear contributions
    // will check this below
    e.setInverted(false);
    // ===========================================
    // non-linear materials
    // ===========================================
    // temporary stress and tangent
    SymmetricMatrix3d sigmaTmp = new SymmetricMatrix3d();
    Matrix6d Dtmp = null;
    if (D != null) {
        Dtmp = new Matrix6d();
    }
    // viscoelastic behaviour
    ViscoelasticBehavior veb = mat.getViscoBehavior();
    double vebTangentScale = 1;
    if (veb != null) {
        vebTangentScale = veb.getTangentScale();
    }
    // incompressibility
    IncompressibleMaterial imat = null;
    if (mat.isIncompressible()) {
        imat = (IncompressibleMaterial) mat;
    }
    MatrixBlock[] constraints = null;
    SymmetricMatrix3d C = new SymmetricMatrix3d();
    // initialize incompressible pressure
    double[] pbuf = myPressures.getBuffer();
    if (mat.isIncompressible() && softIncomp == IncompMethod.ELEMENT) {
        computePressuresAndRinv(e, imat, vebTangentScale);
        if (D != null) {
            constraints = e.getIncompressConstraints();
            for (int i = 0; i < e.myNodes.length; i++) {
                constraints[i].setZero();
            }
        }
    }
    // cache invertible flag
    boolean invertibleMaterials = e.materialsAreInvertible();
    // loop through each integration point
    for (int k = 0; k < ipnts.length; k++) {
        IntegrationPoint3d pt = ipnts[k];
        IntegrationData3d dt = idata[k];
        double scaling = dt.getScaling();
        pt.computeJacobianAndGradient(e.myNodes, idata[k].myInvJ0);
        def.clear();
        def.setF(pt.F);
        def.setAveragePressure(0);
        def.setR(null);
        double detJ = pt.computeInverseJacobian();
        if (detJ < myMinDetJ) {
            myMinDetJ = detJ;
            myMinDetJElement = e;
        }
        if (detJ <= 0 && !invertibleMaterials) {
            e.setInverted(true);
            myNumInverted++;
        }
        // compute shape function gradient and volume fraction
        double dv = detJ * pt.getWeight();
        Vector3d[] GNx = pt.updateShapeGradient(pt.myInvJ);
        // compute pressure
        double pressure = 0;
        double[] H = null;
        if (mat.isIncompressible()) {
            if (softIncomp == IncompMethod.ELEMENT) {
                H = pt.getPressureWeights().getBuffer();
                int npvals = e.numPressureVals();
                for (int l = 0; l < npvals; l++) {
                    pressure += H[l] * pbuf[l];
                }
            } else if (softIncomp == IncompMethod.NODAL) {
                if (e instanceof TetElement) {
                    // use the average pressure for all nodes
                    pressure = 0;
                    for (int i = 0; i < nodes.length; i++) {
                        pressure += nodes[i].myPressure;
                    }
                    pressure /= nodes.length;
                } else if (e.integrationPointsMapToNodes()) {
                    pressure = nodes[k].myPressure;
                } else if (e.integrationPointsInterpolateToNodes()) {
                    // interpolate using shape function
                    VectorNd N = pt.getShapeWeights();
                    pressure = 0;
                    for (int i = 0; i < N.size(); ++i) {
                        pressure += nodes[i].myPressure * N.get(i);
                    }
                }
            } else if (softIncomp == IncompMethod.FULL) {
                pressure = imat.getEffectivePressure(detJ / dt.getDetJ0());
            }
        }
        // anisotropy rotational frame
        Matrix3d Q = (dt.myFrame != null ? dt.myFrame : Matrix3d.IDENTITY);
        // System.out.println("FEM Pressure: " + pressure);
        pt.avgp = pressure;
        def.setAveragePressure(pressure);
        // clear stress/tangents
        pt.sigma.setZero();
        if (D != null) {
            D.setZero();
        }
        // base material
        if (!mat.isLinear()) {
            mat.computeStress(pt.sigma, def, Q, null);
            if (scaling != 1) {
                pt.sigma.scale(scaling);
            }
            if (D != null) {
                mat.computeTangent(D, pt.sigma, def, Q, null);
                if (scaling != 1) {
                    D.scale(scaling);
                }
            }
        }
        // reset pressure to zero for auxiliary materials
        pt.avgp = 0;
        def.setAveragePressure(0);
        if (e.numAuxiliaryMaterials() > 0) {
            for (AuxiliaryMaterial amat : e.getAuxiliaryMaterials()) {
                // skip linear materials
                if (!amat.isLinear()) {
                    amat.computeStress(sigmaTmp, def, pt, dt, mat);
                    if (scaling != 1) {
                        sigmaTmp.scale(scaling);
                    }
                    pt.sigma.add(sigmaTmp);
                    if (D != null) {
                        amat.computeTangent(Dtmp, sigmaTmp, def, pt, dt, mat);
                        if (scaling != 1) {
                            Dtmp.scale(scaling);
                        }
                        D.add(Dtmp);
                    }
                }
            }
        }
        // XXX only uses non-linear stress
        // bring back pressure term
        pt.avgp = pressure;
        def.setAveragePressure(pressure);
        if (veb != null) {
            ViscoelasticState state = idata[k].getViscoState();
            if (state == null) {
                state = veb.createState();
                idata[k].setViscoState(state);
            }
            veb.computeStress(pt.sigma, state);
            if (D != null) {
                veb.computeTangent(D, state);
            }
        } else {
            dt.clearState();
        }
        // sum stress/stiffness contributions to each node
        for (int i = 0; i < e.myNodes.length; i++) {
            FemNode3d nodei = e.myNodes[i];
            int bi = nodei.getSolveIndex();
            FemUtilities.addStressForce(nodei.myInternalForce, GNx[i], pt.sigma, dv);
            if (D != null) {
                double p = 0;
                double kp = 0;
                if (mat.isIncompressible()) {
                    if (softIncomp == IncompMethod.ELEMENT) {
                        FemUtilities.addToIncompressConstraints(constraints[i], H, GNx[i], dv);
                        p = pressure;
                    } else if (softIncomp == IncompMethod.FULL) {
                        double dV = dt.getDetJ0() * pt.getWeight();
                        kp = imat.getEffectiveModulus(detJ / dt.getDetJ0()) * dV;
                        p = pressure;
                    }
                }
                // compute stiffness
                if (bi != -1) {
                    for (int j = 0; j < e.myNodes.length; j++) {
                        int bj = e.myNodes[j].getSolveIndex();
                        if (!mySolveMatrixSymmetricP || bj >= bi) {
                            e.myNbrs[i][j].addMaterialStiffness(GNx[i], D, GNx[j], dv);
                            e.myNbrs[i][j].addGeometricStiffness(GNx[i], pt.sigma, GNx[j], dv);
                            e.myNbrs[i][j].addPressureStiffness(GNx[i], p, GNx[j], dv);
                            if (kp != 0) {
                                e.myNbrs[i][j].addDilationalStiffness(vebTangentScale * kp, GNx[i], GNx[j]);
                            }
                        }
                    }
                }
            }
            // if D != null
            // nodal stress/strain
            double[] nodalExtrapMat = e.getNodalExtrapolationMatrix();
            if (nodalExtrapMat != null) {
                if (myComputeNodalStress) {
                    double a = nodalExtrapMat[i * ipnts.length + k];
                    if (a != 0) {
                        nodei.addScaledStress(a / nodei.numAdjacentElements(), pt.sigma);
                    }
                }
                // if base material non-linear and computing nodal strain
                if (myComputeNodalStrain && !mat.isLinear()) {
                    double a = nodalExtrapMat[i * ipnts.length + k];
                    if (a != 0) {
                        def.computeRightCauchyGreen(C);
                        C.m00 -= 1;
                        C.m11 -= 1;
                        C.m22 -= 1;
                        C.scale(0.5);
                        nodei.addScaledStrain(a / nodei.numAdjacentElements(), C);
                    }
                }
            }
        }
        // nodal incompressibility constraints
        if (D != null && softIncomp == IncompMethod.NODAL && !(e instanceof TetElement)) {
            if (e.integrationPointsMapToNodes()) {
                for (FemNodeNeighbor nbr : getNodeNeighbors(e.myNodes[k])) {
                    int j = e.getLocalNodeIndex(nbr.myNode);
                    if (j != -1) {
                        // myNodalConstraints[i].scale(dv, GNx[i]);
                        nbr.myDivBlk.scaledAdd(dv, GNx[j]);
                    }
                }
            } else if (e.integrationPointsInterpolateToNodes()) {
                // distribute according to shape function weights
                VectorNd N = pt.getShapeWeights();
                for (int i = 0; i < N.size(); ++i) {
                    for (FemNodeNeighbor nbr : getNodeNeighbors(e.myNodes[i])) {
                        int j = e.getLocalNodeIndex(nbr.getNode());
                        if (j != -1) {
                            nbr.myDivBlk.scaledAdd(N.get(i) * dv, GNx[j]);
                        }
                    }
                }
            }
        }
    // soft incompressibility
    }
    // tet nodal incompressibility
    if (D != null && mat.isIncompressible() && softIncomp == IncompMethod.NODAL) {
        if (e instanceof TetElement) {
            ((TetElement) e).getAreaWeightedNormals(myNodalConstraints);
            for (int i = 0; i < 4; i++) {
                myNodalConstraints[i].scale(-1 / 12.0);
            }
            for (int i = 0; i < e.numNodes(); ++i) {
                for (FemNodeNeighbor nbr : getNodeNeighbors(e.myNodes[i])) {
                    int j = e.getLocalNodeIndex(nbr.myNode);
                    if (j != -1) {
                        nbr.myDivBlk.scaledAdd(1, myNodalConstraints[j]);
                    }
                }
            }
        }
    }
    // element-wise incompressibility constraints
    if (D != null) {
        if (mat.isIncompressible() && softIncomp == IncompMethod.ELEMENT) {
            boolean kpIsNonzero = false;
            int npvals = e.numPressureVals();
            for (int l = 0; l < npvals; l++) {
                double Jpartial = e.myVolumes[l] / e.myRestVolumes[l];
                myKp[l] = imat.getEffectiveModulus(Jpartial) / e.myRestVolumes[l];
                if (myKp[l] != 0) {
                    kpIsNonzero = true;
                }
            }
            // double kp = imat.getEffectiveModulus(vol/restVol)/restVol;
            if (true) {
                for (int i = 0; i < e.myNodes.length; i++) {
                    int bi = e.myNodes[i].getSolveIndex();
                    if (bi != -1) {
                        for (int j = 0; j < e.myNodes.length; j++) {
                            int bj = e.myNodes[j].getSolveIndex();
                            if (!mySolveMatrixSymmetricP || bj >= bi) {
                                e.myNbrs[i][j].addDilationalStiffness(myRinv, constraints[i], constraints[j]);
                            }
                        // end filling in symmetric
                        }
                    // end filling in dilatational stiffness
                    }
                // end checking if valid index
                }
            // end looping through nodes
            }
        // XXX ALWAYS??
        }
    // end soft elem incompress
    }
// end checking if computing tangent
}

Example 25 with SymmetricMatrix3d

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

the class FemNode3d method setStrain.

public void setStrain(double vms) {
    if (myAvgStrain == null) {
        myAvgStrain = new SymmetricMatrix3d();
    }
    myAvgStrain.setZero();
    myAvgStrain.m00 = vms;
    myAvgStrain.m11 = vms;
    myAvgStrain.m22 = vms;
}