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

Example 6 with FemMaterial

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

the class AuxMaterialElementDesc method computeTangent.

@Override
public void computeTangent(Matrix6d D, SymmetricMatrix3d stress, 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;
            if (Q == null) {
                Q = Matrix3d.IDENTITY;
            }
            mat.computeTangent(D, stress, def, Q, baseMat);
            D.scale(frac);
        } else {
            D.setZero();
        }
    }
}

Example 7 with FemMaterial

use of artisynth.core.materials.FemMaterial 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 8 with FemMaterial

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

the class FemMuscleHeart method build.

// Model builder
@Override
public void build(String[] args) throws IOException {
    super.build(args);
    setMaxStepSize(0.005);
    // Root mechanical model
    MechModel mech = new MechModel("mech");
    mech.setGravity(0, 0, -9.8);
    addModel(mech);
    // -------------------------------------------------------------
    // HEART LOAD / ADD GEOMETRY
    // -------------------------------------------------------------
    // Heart surface mesh, with texture
    String heartFile = ArtisynthPath.getSrcRelativePath(this, "data/HumanHeart.obj");
    WavefrontReader wfr = new WavefrontReader(new File(heartFile));
    PolygonalMesh heartMesh = new PolygonalMesh();
    wfr.readMesh(heartMesh);
    // triangulate for interaction
    heartMesh.triangulate();
    // FEM heart:
    // - FEM mesh of heart convex hull
    // - embedded heart surface geometry
    FemMuscleModel heart = new FemMuscleModel("heart");
    TetGenReader.read(heart, ArtisynthPath.getSrcRelativePath(this, "data/HumanHeartHull.node"), ArtisynthPath.getSrcRelativePath(this, "data/HumanHeartHull.ele"));
    // add real-looking mesh
    FemMeshComp embeddedHeart = heart.addMesh(heartMesh);
    embeddedHeart.setName("embedded");
    // Allow inverted elements (poor quality mesh)
    heart.setWarnOnInvertedElements(false);
    heart.setAbortOnInvertedElements(false);
    // Convert unites to metres (original was cm)
    heart.scaleDistance(0.01);
    heart.setGravity(0, 0, -9.8);
    heart.setStiffnessDamping(0.02);
    // Set material properties
    heart.setDensity(1000);
    FemMaterial femMat = new LinearMaterial(2500, 0.33, true);
    // simple muscle
    MuscleMaterial muscleMat = new SimpleForceMuscle(500.0);
    heart.setMaterial(femMat);
    // Add heart to model
    mech.addModel(heart);
    // -------------------------------------------------------------
    // MUSCLE BUNDLES
    // -------------------------------------------------------------
    // One "long" direction muscle bundle
    // One "radial" muscle bundle
    // LONG BUNDLE
    // Compute the "long" direction of the heart
    PolygonalMesh hull = heart.getSurfaceMesh();
    RigidTransform3d trans = hull.computePrincipalAxes();
    Vector3d longAxis = new Vector3d();
    // first column of rotation
    trans.R.getColumn(0, longAxis);
    // Create the long axis muscle bundle
    MuscleBundle longBundle = new MuscleBundle("long");
    for (FemElement3d elem : heart.getElements()) {
        longBundle.addElement(elem, longAxis);
    }
    longBundle.setMuscleMaterial(muscleMat);
    heart.addMuscleBundle(longBundle);
    // RADIAL BUNDLE
    // Compute a plane through centre of heart
    Plane plane = new Plane(longAxis, new Point3d(trans.p));
    Point3d centroid = new Point3d();
    Vector3d radialDir = new Vector3d();
    // Create the radial muscle bundle
    MuscleBundle radialBundle = new MuscleBundle("radial");
    for (FemElement3d elem : heart.getElements()) {
        elem.computeCentroid(centroid);
        // project to plane and compute radial direction
        plane.project(centroid, centroid);
        radialDir.sub(centroid, trans.p);
        radialDir.normalize();
        radialBundle.addElement(elem, radialDir);
    }
    radialBundle.setMuscleMaterial(muscleMat);
    heart.addMuscleBundle(radialBundle);
    // -------------------------------------------------------------
    // RIGID TABLE AND COLLISION
    // -------------------------------------------------------------
    // Create a rigid box for the heart to fall on
    RigidBody box = RigidBody.createBox("box", 0.2, 0.2, 0.02, 0, /*addnormals*/
    true);
    box.setPose(new RigidTransform3d(new Vector3d(0, 0, -0.2), AxisAngle.IDENTITY));
    box.setDynamic(false);
    mech.addRigidBody(box);
    // Enable collisions between the heart and table
    mech.setCollisionBehavior(heart, box, true);
    // -------------------------------------------------------------
    // RENDER PROPERTIES
    // -------------------------------------------------------------
    // Hide elements and nodes
    RenderProps.setVisible(heart.getElements(), false);
    RenderProps.setVisible(heart.getNodes(), false);
    RenderProps.setLineColor(radialBundle, Color.BLUE);
    RenderProps.setLineColor(longBundle, Color.RED);
    radialBundle.setDirectionRenderLen(0.1);
    longBundle.setDirectionRenderLen(0.1);
    RenderProps.setVisible(radialBundle, false);
    RenderProps.setVisible(longBundle, false);
    RenderProps.setVisible(heart.getSurfaceMeshComp(), false);
    // adjust table render properties
    RenderProps.setShading(box, Shading.METAL);
    RenderProps.setSpecular(box, new Color(0.8f, 0.8f, 0.8f));
    // adjust heart mesh render properties
    RenderProps rprops = embeddedHeart.getRenderProps();
    rprops.getBumpMap().setScaling(0.01f);
    // don't modify specular
    rprops.getColorMap().setSpecularColoring(false);
    rprops.setShading(Shading.SMOOTH);
    rprops.setFaceColor(new Color(0.8f, 0.8f, 0.8f));
    rprops.getColorMap().setColorMixing(ColorMixing.MODULATE);
    rprops.setSpecular(new Color(0.4f, 0.4f, 0.4f));
    rprops.setShininess(128);
    // -------------------------------------------------------------
    // INPUT PROBES
    // -------------------------------------------------------------
    // Add heart probe
    addHeartProbe(longBundle, radialBundle);
}

Example 9 with FemMaterial

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

the class FemModel3d method updateStressAndStiffness.

// DIVBLK
public void updateStressAndStiffness() {
    // allocate or deallocate nodal incompressibility blocks
    setNodalIncompBlocksAllocated(getSoftIncompMethod() == IncompMethod.NODAL);
    // clear existing internal forces and maybe stiffnesses
    for (FemNode3d n : myNodes) {
        n.myInternalForce.setZero();
        if (!myStiffnessesValidP) {
            for (FemNodeNeighbor nbr : getNodeNeighbors(n)) {
                nbr.zeroStiffness();
            }
            // used for soft nodal-based incompressibilty:
            for (FemNodeNeighbor nbr : getIndirectNeighbors(n)) {
                nbr.zeroStiffness();
            }
        }
        if (myComputeNodalStress) {
            n.zeroStress();
        }
        if (myComputeNodalStrain) {
            n.zeroStrain();
        }
    }
    if (!myVolumeValid) {
        updateVolume();
    }
    IncompMethod softIncomp = getSoftIncompMethod();
    if (softIncomp == IncompMethod.NODAL) {
        if (!myNodalRestVolumesValidP) {
            updateNodalRestVolumes();
        }
        setNodalIncompConstraintsAllocated(true);
        updateNodalPressures((IncompressibleMaterial) myMaterial);
        for (FemNode3d n : myNodes) {
            for (FemNodeNeighbor nbr : getNodeNeighbors(n)) {
                nbr.myDivBlk.setZero();
            }
        }
    }
    Matrix6d D = new Matrix6d();
    // compute new forces as well as stiffness matrix if warping is enabled
    myMinDetJ = Double.MAX_VALUE;
    myMinDetJElement = null;
    myNumInverted = 0;
    double mins = Double.MAX_VALUE;
    FemElement3d minE = null;
    for (FemElement3d e : myElements) {
        FemMaterial mat = getElementMaterial(e);
        computeStressAndStiffness(e, mat, D, softIncomp);
        if (checkTangentStability) {
            double s = checkMatrixStability(D);
            if (s < mins) {
                mins = s;
                minE = e;
            }
        }
    }
    // incompressibility
    if ((softIncomp == IncompMethod.NODAL) && myMaterial != null && myMaterial.isIncompressible()) {
        computeNodalIncompressibility((IncompressibleMaterial) myMaterial, D);
    }
    if (checkTangentStability && minE != null) {
        System.out.println("min s=" + mins + ", element " + minE.getNumber());
    }
    if (myNumInverted > 0) {
        if (myWarnOnInvertedElems) {
            System.out.println("Warning: " + myNumInverted + " inverted elements; min detJ=" + myMinDetJ + ", element " + ComponentUtils.getPathName(myMinDetJElement));
        }
        if (myAbortOnInvertedElems) {
            throw new NumericalException("Inverted elements");
        }
    }
    if (!myStiffnessesValidP && mySolveMatrixSymmetricP) {
        for (FemNode3d n : myNodes) {
            int bi = n.getSolveIndex();
            if (bi != -1) {
                for (FemNodeNeighbor nbr : getNodeNeighbors(n)) {
                    int bj = nbr.myNode.getSolveIndex();
                    if (bj > bi) {
                        FemNodeNeighbor nbrT = nbr.myNode.getNodeNeighborBySolveIndex(bi);
                        nbrT.setTransposedStiffness(nbr);
                    }
                }
                // used for soft nodal-based incompressibilty:
                for (FemNodeNeighbor nbr : getIndirectNeighbors(n)) {
                    int bj = nbr.myNode.getSolveIndex();
                    if (bj > bi) {
                        FemNodeNeighbor nbrT = nbr.myNode.getIndirectNeighborBySolveIndex(bi);
                        nbrT.setTransposedStiffness(nbr);
                    }
                }
            }
        }
    }
    myStiffnessesValidP = true;
    myStressesValidP = true;
// timerStop("stressAndStiffness");
}

Example 10 with FemMaterial

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

the class FemModel3d method getNodalDeformationGradients.

public void getNodalDeformationGradients(Matrix3d[] Fnodal) {
    if (Fnodal.length < myNodes.size()) {
        throw new IllegalArgumentException("Fnodal must have length >= " + myNodes.size());
    }
    for (FemElement3d e : myElements) {
        FemNode3d[] enodes = e.myNodes;
        FemMaterial mat = getElementMaterial(e);
        if (mat.isLinear()) {
            IntegrationPoint3d wpnt = e.getWarpingPoint();
            IntegrationData3d data = e.getWarpingData();
            wpnt.computeJacobianAndGradient(enodes, data.myInvJ0);
            for (int i = 0; i < enodes.length; i++) {
                int nidx = myNodes.indexOf(enodes[i]);
                Fnodal[nidx].scaledAdd(1.0 / enodes[i].numAdjacentElements(), wpnt.F);
            }
        } else {
            IntegrationPoint3d[] ipnts = e.getIntegrationPoints();
            IntegrationData3d[] idata = e.getIntegrationData();
            double[] nodalExtrapMat = e.getNodalExtrapolationMatrix();
            for (int k = 0; k < ipnts.length; k++) {
                ipnts[k].computeJacobianAndGradient(e.myNodes, idata[k].myInvJ0);
                for (int i = 0; i < enodes.length; i++) {
                    double a = nodalExtrapMat[i * ipnts.length + k];
                    int nidx = myNodes.indexOf(enodes[i]);
                    Fnodal[nidx].scaledAdd(a / enodes[i].numAdjacentElements(), ipnts[k].F);
                }
            }
        }
    }
}