Examples with Point artisynth.core.mechmodels.Point used on opensource projects

Example 21 with Point

use of artisynth.core.mechmodels.Point in project artisynth_core by artisynth.

the class FemElement3d method isInvertedAtRest.

/**
 * Returns true if the rest position for this element results in a negative
 * Jacobian determinant for at least one integration point. This method
 * should be used sparingly since it is computes this result from scratch.
 *
 * @return true if the rest position is inverted.
 */
public boolean isInvertedAtRest() {
    IntegrationPoint3d[] ipnts = getIntegrationPoints();
    for (int i = 0; i < ipnts.length; i++) {
        IntegrationData3d idata = new IntegrationData3d();
        /*
          *  code snippet from IntegrationData3d.computeRestJacobian()
          *  -- there may be better way to check "invertedness"
          */
        Matrix3d J0 = new Matrix3d();
        J0.setZero();
        for (int j = 0; j < myNodes.length; j++) {
            Vector3d pos = myNodes[j].myRest;
            Vector3d dNds = ipnts[i].GNs[j];
            J0.addOuterProduct(pos.x, pos.y, pos.z, dNds.x, dNds.y, dNds.z);
        }
        if (idata.getInvJ0().fastInvert(J0) <= 0) {
            return true;
        }
    }
    return false;
}

Example 22 with Point

use of artisynth.core.mechmodels.Point in project artisynth_core by artisynth.

the class TextLabeller3d method labelPoints.

public void labelPoints(List<? extends Point> pnts) {
    for (Point pnt : pnts) {
        String text = pnt.getName();
        if (text == null) {
            text = "" + pnt.getNumber();
        }
        addItem(text, pnt.getPosition(), true);
    }
}

Example 23 with Point

use of artisynth.core.mechmodels.Point in project artisynth_core by artisynth.

the class TextLabeller3d method labelPoints.

public void labelPoints(PointList<? extends Point> pnts) {
    for (Point pnt : pnts) {
        String text = pnt.getName();
        if (text == null) {
            text = "" + pnt.getNumber();
        }
        addItem(text, pnt.getPosition(), true);
    }
}

Example 24 with Point

use of artisynth.core.mechmodels.Point 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 Point

use of artisynth.core.mechmodels.Point in project artisynth_core by artisynth.

the class FemMeshComp method createEmbedded.

public static FemMeshComp createEmbedded(FemMeshComp surf, MeshBase mesh, FemModel3d fem) {
    double reduceTol = 1e-8;
    ArrayList<FemNode> nodes = new ArrayList<FemNode>();
    VectorNd weights = new VectorNd();
    if (surf == null) {
        surf = new FemMeshComp(fem);
    }
    surf.setMesh(mesh);
    ArrayList<Vertex3d> verts = mesh.getVertices();
    surf.myVertexAttachments.clear();
    for (int i = 0; i < verts.size(); i++) {
        // this could works very similarly to the code that adds
        // marker points into a mesh
        Vertex3d vtx = verts.get(i);
        // if (vtx instanceof FemMeshVertex) {
        // nodes.clear();
        // nodes.add(((FemMeshVertex)vtx).getPoint());
        // weights.clear();
        // weights.add(1.0);
        // }
        // else
        {
            FemElement3d elem = surf.myFem.findContainingElement(vtx.pnt);
            Point3d newLoc = new Point3d(vtx.pnt);
            if (elem == null) {
                // won't use newLoc since we're not projecting vertex onto FEM
                elem = surf.myFem.findNearestSurfaceElement(newLoc, vtx.pnt);
            }
            VectorNd coords = new VectorNd(elem.numNodes());
            // first see if there's a node within reduceTol of the point,
            // and if so just use that
            double maxDist = Double.NEGATIVE_INFINITY;
            double minDist = Double.POSITIVE_INFINITY;
            FemNode3d nearestNode = null;
            FemNode3d[] elemNodes = elem.getNodes();
            for (int k = 0; k < elemNodes.length; k++) {
                double d = vtx.pnt.distance(elemNodes[k].getPosition());
                if (d > maxDist) {
                    maxDist = d;
                }
                if (d < minDist) {
                    minDist = d;
                    nearestNode = elemNodes[k];
                }
            }
            if (minDist / maxDist <= reduceTol) {
                // weight everything to the nearest node
                nodes.clear();
                nodes.add(nearestNode);
                weights.setSize(0);
                weights.append(1.0);
            } else {
                Vector3d c3 = new Vector3d();
                boolean converged = elem.getNaturalCoordinates(c3, vtx.pnt, 1000) >= 0;
                if (!converged) {
                    System.err.println("Warning: getNaturalCoordinatesRobust() did not converge, " + "element=" + ComponentUtils.getPathName(elem) + ", point=" + vtx.pnt);
                // c3.setZero();
                // XXX debugging:
                // elem.getNaturalCoordinates(c3,  vtx.pnt, 1000); // try again once more
                }
                for (int j = 0; j < elem.numNodes(); j++) {
                    coords.set(j, elem.getN(j, c3));
                }
                nodes.clear();
                weights.setSize(0);
                for (int k = 0; k < coords.size(); k++) {
                    if (Math.abs(coords.get(k)) >= reduceTol) {
                        nodes.add(elem.getNodes()[k]);
                        weights.append(coords.get(k));
                    }
                }
            }
        }
        if (weights.size() > 1) {
            PointFem3dAttachment attacher = new PointFem3dAttachment();
            attacher.setFromNodes(nodes, weights);
            surf.myVertexAttachments.add(attacher);
        } else if (weights.size() == 1) {
            PointParticleAttachment attacher = new PointParticleAttachment(nodes.get(0), null);
            surf.myVertexAttachments.add(attacher);
        }
    }
    surf.buildNodeVertexMap();
    return surf;
}