Examples with ViscoelasticState artisynth.core.materials.ViscoelasticState used on opensource projects

Example 1 with ViscoelasticState

use of artisynth.core.materials.ViscoelasticState in project artisynth_core by artisynth.

the class FemModel3d method advanceAuxState.

public void advanceAuxState(double t0, double t1) {
    for (int i = 0; i < myElements.size(); i++) {
        FemElement3d e = myElements.get(i);
        FemMaterial mat = getElementMaterial(e);
        if (mat.getViscoBehavior() != null) {
            ViscoelasticBehavior veb = mat.getViscoBehavior();
            IntegrationData3d[] idata = e.getIntegrationData();
            for (int k = 0; k < idata.length; k++) {
                ViscoelasticState state = idata[k].getViscoState();
                if (state == null) {
                    state = veb.createState();
                    idata[k].setViscoState(state);
                }
                veb.advanceState(state, t0, t1);
            }
        }
    }
}

Example 2 with ViscoelasticState

use of artisynth.core.materials.ViscoelasticState 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
}