Tutorial: Non-Linear solid analysis using the gsStructuralAnalysis module

Static simulations of a solid

Annotated source file

Here is the full file nonlinear_solid_static.cpp. Clicking on a function or class name will lead you to its reference documentation.

 
#include <gismo.h>
 
#ifdef gsElasticity_ENABLED
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
#include <gsElasticity/gsGeoUtils.h>
#endif
 
 
#include <gsStructuralAnalysis/src/gsStaticSolvers/gsStaticNewton.h>
#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisUtils.h>
 
using namespace gismo;
 
#ifdef gsElasticity_ENABLED
int main(int argc, char *argv[])
{
    bool plot  = false;
    index_t numRefine  = 1;
    index_t numElevate = 1;
 
    gsCmdLine cmd("Linear shell tutorial.");
    cmd.addInt( "e", "degreeElevation","Number of degree elevation steps to perform before solving (0: equalize degree in all directions)", numElevate );
    cmd.addInt( "r", "uniformRefine", "Number of Uniform h-refinement steps to perform before solving",  numRefine );
    cmd.addSwitch("plot", "Create a ParaView visualization file with the solution", plot);
    try { cmd.getValues(argc,argv); } catch (int rv) { return rv; }
 
    // Initialize [ori]ginal and [def]ormed geometry
    gsMultiPatch<> ori;
    std::string fileName = "paraboloid_volume.xml";
    gsReadFile<>(fileName, ori);
 
    // p-refine
    if (numElevate!=0)
        ori.degreeElevate(numElevate);
    // h-refine (only in first two directions)
    for (int r =0; r < numRefine; ++r)
    {
        ori.uniformRefine(1,1,0);
        ori.uniformRefine(1,1,1);
    }
 
    // creating basis
    gsMultiBasis<> basis(ori);
 
    // Define the boundary conditions object
    gsBoundaryConditions<> bc;
    // Set the geometry map for computation of the Dirichlet BCs
    bc.setGeoMap(ori);
 
    // Set the boundary conditions
    gsFunctionExpr<> surf_force("0","0","-1e5",3);
    for (index_t c=0; c!=3; c++)
    {
        bc.addCornerValue(boundary::southwestfront, 0.0, 0, c); // (corner,value, patch, unknown)
        bc.addCornerValue(boundary::southeastfront, 0.0, 0, c); // (corner,value, patch, unknown)
        bc.addCornerValue(boundary::northwestfront, 0.0, 0, c); // (corner,value, patch, unknown)
        bc.addCornerValue(boundary::northeastfront, 0.0, 0, c); // (corner,value, patch, unknown)
    }
    bc.addCondition(boundary::front, condition_type::neumann, &surf_force);
    bc.addCondition(boundary::front, condition_type::neumann, &surf_force);
 
    // The surface force is defined in the physical space, i.e. 3D
    gsFunctionExpr<> body_force("0","0","0",3);
 
    gsElasticityAssembler<real_t> assembler(ori,basis,bc,body_force);
    assembler.options().setReal("YoungsModulus",1.E9);
    assembler.options().setReal("PoissonsRatio",0.45);
    assembler.options().setInt("MaterialLaw",material_law::saint_venant_kirchhoff);
    assembler.options().setInt("DirichletValues",dirichlet::l2Projection);
 
    std::vector<gsMatrix<> > fixedDofs = assembler.allFixedDofs();
    // Function for the Jacobian
    gsStructuralAnalysisOps<real_t>::Jacobian_t Jacobian = [&assembler,&fixedDofs](gsVector<real_t> const &x, gsSparseMatrix<real_t> & m)
    {
        assembler.assemble(x,fixedDofs);
        m = assembler.matrix();
        return true;
    };
 
    // Function for the Residual
    gsStructuralAnalysisOps<real_t>::Residual_t Residual = [&fixedDofs,&assembler](gsVector<real_t> const &x, gsVector<real_t> & v)
    {
        assembler.assemble(x,fixedDofs);
        v = assembler.rhs();
        return true;
    };
 
    assembler.assemble();
    gsSparseMatrix<> K = assembler.matrix();
    gsVector<> F = assembler.rhs();
 
    gsStaticNewton<real_t> solver(K,F,Jacobian,Residual);
    solver.options().setInt("verbose",1);
    solver.initialize();
 
    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsStatus status = solver.solve();
    GISMO_ASSERT(status==gsStatus::Success,"Solver failed");
    gsVector<> solVector = solver.solution();
 
    // constructing solution as an IGA function
    gsMultiPatch<> displ;
    assembler.constructSolution(solVector,assembler.allFixedDofs(),displ);
    gsMultiPatch<> def = ori;
    for (index_t p = 0; p < def.nPieces(); ++p)
        def.patch(p).coefs() += displ.patch(p).coefs();
 
    // Evaluate the solution on a reference point (parametric coordinates)
    // The reference points are stored column-wise
    gsMatrix<> refPoint(3,1);
    refPoint<<0.5,0.5,1.0;
    // Compute the values
    gsVector<> physpoint = def.patch(0).eval(refPoint);
    gsVector<> refVals = displ.patch(0).eval(refPoint);
    // gsInfo << "Displacement at reference point: "<<numVal<<"\n";
    gsInfo  << "Displacement at reference point ("
            <<ori.patch(0).eval(refPoint).transpose()<<"): "
            <<": "<<refVals.transpose()<<"\n";
 
    // ! [Export visualization in ParaView]
    if (plot)
    {
        // Plot the displacements on the deformed geometry
        gsField<> solField(def, displ);
        gsInfo<<"Plotting in Paraview...\n";
        gsWriteParaview<>( solField, "Deformation", 10000, true);
 
        // Plot the membrane Von-Mises stresses on the geometry
        gsPiecewiseFunction<> VMm;
        assembler.constructCauchyStresses(displ,VMm,stress_components::von_mises);
        gsField<> membraneStress(def,VMm, true);
        gsWriteParaview(membraneStress,"MembraneStress",10000);
    }
    // ! [Export visualization in ParaView]
 
    return EXIT_SUCCESS;
#else
    GISMO_ERROR("The tutorial needs to be compiled with gsElasticity enabled");
    return EXIT_FAILED;
#endif
 
}// end main