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

Example 1 with Matrix

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

the class FemModel3d method computePressuresAndRinv.

protected void computePressuresAndRinv(FemElement3d e, IncompressibleMaterial imat, double scale) {
    int npvals = e.numPressureVals();
    myRinv.setSize(npvals, npvals);
    myPressures.setSize(npvals);
    double[] pbuf = myPressures.getBuffer();
    double restVol = e.getRestVolume();
    if (npvals > 1) {
        myPressures.setZero();
        IntegrationPoint3d[] ipnts = e.getIntegrationPoints();
        IntegrationData3d[] idata = e.getIntegrationData();
        if (imat.getBulkPotential() != BulkPotential.QUADRATIC) {
            myRinv.setZero();
        }
        for (int k = 0; k < ipnts.length; k++) {
            IntegrationPoint3d pt = ipnts[k];
            pt.computeJacobian(e.getNodes());
            double detJ0 = idata[k].getDetJ0();
            double detJ = pt.getJ().determinant() / detJ0;
            double dV = detJ0 * pt.getWeight();
            double[] H = pt.getPressureWeights().getBuffer();
            for (int i = 0; i < npvals; i++) {
                pbuf[i] += H[i] * imat.getEffectivePressure(detJ) * dV;
            }
            if (imat.getBulkPotential() != BulkPotential.QUADRATIC) {
                double mod = imat.getEffectiveModulus(detJ);
                for (int i = 0; i < npvals; i++) {
                    for (int j = 0; j < npvals; j++) {
                        myRinv.add(i, j, H[i] * H[j] * mod * dV);
                    }
                }
            }
        }
        Matrix W = e.getPressureWeightMatrix();
        W.mul(myPressures, myPressures);
        myPressures.scale(1 / restVol);
        if (imat.getBulkPotential() == BulkPotential.QUADRATIC) {
            myRinv.set(W);
            myRinv.scale(scale * imat.getBulkModulus() / restVol);
        } else {
            // optimize later
            MatrixNd Wtmp = new MatrixNd(W);
            Wtmp.scale(scale / restVol);
            myRinv.mul(Wtmp);
            myRinv.mul(Wtmp, myRinv);
        }
    } else {
        double Jpartial = e.myVolumes[0] / e.myRestVolumes[0];
        pbuf[0] = (imat.getEffectivePressure(Jpartial) + 0 * e.myLagrangePressures[0]);
        myRinv.set(0, 0, scale * imat.getEffectiveModulus(Jpartial) / restVol);
    }
}