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

Example 11 with MatrixBlock

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

the class DynamicAttachment method reduceRowMatrix.

/**
 * Reduces a sparse row matrix (such as a constraint
 * matrix) to account for this attachment. This is done by applying
 * the transform
 * <pre>
 *        T
 * G = G P
 * </pre>
 * where
 * <pre>
 *             T
 *     [ I -Gka  ]
 * P = [         ]
 *     [ 0   0   ]
 * </pre>
 * and Gka is the constraint matrix for this attachment.
 *
 * <p>At present, this method requires the matrix G to be vertically linked.
 */
public void reduceRowMatrix(SparseBlockMatrix G) {
    if (!G.isVerticallyLinked()) {
        throw new IllegalArgumentException("Matrix G is not vertically linked");
    }
    DynamicComponent[] masters = getMasters();
    if (masters.length == 0) {
        return;
    }
    int bs = getSlaveSolveIndex();
    if (bs == -1) {
        return;
    }
    MatrixBlock blk = G.firstBlockInCol(bs);
    while (blk != null) {
        int bi = blk.getBlockRow();
        for (int j = 0; j < masters.length; j++) {
            if (masters[j].getSolveIndex() != -1) {
                int bm = masters[j].getSolveIndex();
                MatrixBlock depBlk = G.getBlock(bi, bm);
                if (depBlk == null) {
                    depBlk = MatrixBlockBase.alloc(blk.rowSize(), G.getBlockColSize(bm));
                    // depBlk = createColBlock (blk.rowSize());
                    G.addBlock(bi, bm, depBlk);
                }
                mulSubG(depBlk, blk, j);
            }
        }
        blk = blk.down();
    }
}

Example 12 with MatrixBlock

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

the class DynamicAttachment method addAttachmentJacobian.

/**
 * Reduces the system matrix to account for this attachment. This
 * is done by applying the transform
 * <pre>
 *          T
 * M = P M P
 * </pre>
 * where
 * <pre>
 *             T
 *     [ I -Gka  ]
 * P = [         ]
 *     [ 0   0   ]
 * </pre>
 * and Gka is the constraint matrix for this attachment.
 */
public void addAttachmentJacobian(SparseBlockMatrix S, VectorNd f, boolean[] reduced) {
    DynamicComponent[] masters = getMasters();
    if (masters.length == 0) {
        return;
    }
    if (myMasterBlks == null || masters.length > myMasterBlks.length) {
        // System.out.println ("new num master blocks: " + masters.length);
        myMasterBlks = new MatrixBlock[masters.length];
    }
    int bm;
    int bs = getSlaveSolveIndex();
    // System.out.println ("bs=" + bs);
    if (bs == -1) {
        return;
    }
    // boolean debug =
    // (this instanceof PointParticleAttachment && getSlave().getNumber() == 0);
    double[] dbuf = null;
    double[] fbuf = null;
    int soff = -1;
    int ssize = 0;
    if (f != null) {
        fbuf = f.getBuffer();
        soff = S.getBlockRowOffset(bs);
        ssize = S.getBlockRowSize(bs);
        dbuf = getNegatedDerivative(new VectorNd(ssize));
    }
    // 
    // rows first: M = P M
    // 
    MatrixBlock blk = S.firstBlockInRow(bs);
    // get first row block for each master
    for (int i = 0; i < masters.length; i++) {
        if ((bm = masters[i].getSolveIndex()) != -1) {
            myMasterBlks[i] = S.firstBlockInRow(bm);
        }
    }
    while (blk != null) {
        int bj = blk.getBlockCol();
        for (int i = 0; i < masters.length; i++) {
            if ((bm = masters[i].getSolveIndex()) != -1) {
                MatrixBlock depBlk = myMasterBlks[i];
                while (depBlk.getBlockCol() < bj) {
                    depBlk = depBlk.next();
                }
                if (depBlk.getBlockCol() != bj) {
                    throw new InternalErrorException("slave blk at (" + bs + "," + bj + "), master at (" + depBlk.getBlockRow() + "," + depBlk.getBlockCol() + ")");
                }
                mulSubGT(depBlk, blk, i);
                myMasterBlks[i] = depBlk;
            }
        }
        blk.setZero();
        blk = blk.next();
    }
    reduced[bs] = true;
    // T
    // columns next: M = M P, and f += M dg
    // 
    blk = S.firstBlockInCol(bs);
    for (int i = 0; i < masters.length; i++) {
        if ((bm = masters[i].getSolveIndex()) != -1) {
            myMasterBlks[i] = S.firstBlockInCol(bm);
        }
    }
    while (blk != null) {
        int bi = blk.getBlockRow();
        if (!reduced[bi]) {
            for (int i = 0; i < masters.length; i++) {
                if ((bm = masters[i].getSolveIndex()) != -1) {
                    MatrixBlock depBlk = myMasterBlks[i];
                    while (depBlk.getBlockRow() < bi) {
                        depBlk = depBlk.down();
                    }
                    if (depBlk.getBlockRow() != bi) {
                        throw new InternalErrorException("slave blk at (" + bi + "," + bs + "), master at (" + depBlk.getBlockRow() + "," + depBlk.getBlockCol() + ")");
                    }
                    mulSubG(depBlk, blk, i);
                    myMasterBlks[i] = depBlk;
                }
            }
            if (fbuf != null && dbuf != null) {
                int foff = S.getBlockRowOffset(bi);
                blk.mulAdd(fbuf, foff, dbuf, 0);
            }
            blk.setZero();
        }
        blk = blk.down();
    }
    if (fbuf != null) {
        // f = P f
        for (int i = 0; i < masters.length; i++) {
            if ((bm = masters[i].getSolveIndex()) != -1) {
                int moff = S.getBlockRowOffset(bm);
                mulSubGT(fbuf, moff, fbuf, soff, i);
            }
        }
        // zero out the slave term
        for (int i = 0; i < ssize; i++) {
            fbuf[soff + i] = 0;
        }
    }
}

Example 13 with MatrixBlock

use of maspack.matrix.MatrixBlock 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 14 with MatrixBlock

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

the class FemModel3d method computeAvgGNx.

/**
 * Computes the average deformation gradient for an element.
 */
protected void computeAvgGNx(FemElement3d e) {
    IntegrationPoint3d[] ipnts = e.getIntegrationPoints();
    IntegrationData3d[] idata = e.getIntegrationData();
    Vector3d[] avgGNx = null;
    MatrixBlock[] constraints = null;
    constraints = e.getIncompressConstraints();
    for (int i = 0; i < e.myNodes.length; i++) {
        constraints[i].setZero();
    }
    e.setInverted(false);
    for (int k = 0; k < ipnts.length; k++) {
        IntegrationPoint3d pt = ipnts[k];
        pt.computeJacobianAndGradient(e.myNodes, idata[k].myInvJ0);
        double detJ = pt.computeInverseJacobian();
        if (detJ <= 0) {
            e.setInverted(true);
        // if (abortOnInvertedElems) {
        // throw new NumericalException ("Inverted elements");
        // }
        }
        double dv = detJ * pt.getWeight();
        Vector3d[] GNx = pt.updateShapeGradient(pt.myInvJ);
        double[] H = pt.getPressureWeights().getBuffer();
        for (int i = 0; i < e.myNodes.length; i++) {
            FemUtilities.addToIncompressConstraints(constraints[i], H, GNx[i], dv);
        }
    }
}

Example 15 with MatrixBlock

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

the class PointFem3dAttachment method mulSubG.

public void mulSubG(MatrixBlock D, MatrixBlock B, int idx) {
    if (myNoFrameRelativeP) {
        addBlock(D, B, idx);
    } else {
        MatrixBlock G = myMasterBlocks[idx].createTranspose();
        D.mulAdd(B, G);
    }
}