Tutorial: Non-Linear Kirchhoff-Love shell analysis using the gsStructuralAnalysis module

This example solves a non-linear shell problem using the gsThinShellAssembler and a static solver from Structural Analysis Module. The goals of the tutorial are the following:

• Solve a geometrically non-linear problem using gsStaticNewton

For detailed information about the construction of gsThinShellAssembler and gsMaterialMatrixBase, please consult the tutorials Tutorial: Linear Kirchhoff-Love shell analysis and Tutorial: Non-Linear Kirchhoff-Love shell analysis.

Problem set-up

The set-up of the problem to be solved for this example is identical to the one in Tutorial: Non-Linear Kirchhoff-Love shell analysis. The reader is referred to that example for a full overview of the code.

Jacobian and Residual

In Tutorial: Non-Linear Kirchhoff-Love shell analysis, the Jacobian and Residual operators are already defined as std::function objects. In the present example, we use a similar approach, now defined through gsStructuralAnalysisOps. In particular, the Jacobian and Residual operators are given by:

    // Function for the Jacobian
    gsStructuralAnalysisOps<real_t>::Jacobian_t Jacobian = [&assembler](gsVector<real_t> const &x, gsSparseMatrix<real_t> & m)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleMatrix(def);
        m = assembler.matrix();
        return status==ThinShellAssemblerStatus::Success;
    };
    // Function for the Residual
    gsStructuralAnalysisOps<real_t>::Residual_t Residual = [&assembler](gsVector<real_t> const &x, gsVector<real_t> & v)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleVector(def);
        v = assembler.rhs();
        return status==ThinShellAssemblerStatus::Success;
    };

Define the static solver

The Newton-Raphson solver - gsStaticNewton - used in this example is derived from gsStaticBase. There are also other solvers available in the Static solvers submodule. The Newton-Raphson solver needs to be defined using a linear stiffness matrix \(K\), an external force vector \(F\), a Jacobian matrix \(K_{NL}(\mathbf{u})\) and a residual vector \(R(\mathbf{u})\). The latter two are defined above in Jacobian and Residual, whereas the former two are defined like in linear_shell_Assembler:

    assembler.assemble();
    gsSparseMatrix<> K = assembler.matrix();
    gsVector<> F = assembler.rhs();

Consequently, the static solver is constructed as follows:

    gsStaticNewton<real_t> solver(K,F,Jacobian,Residual);
    solver.options().setInt("verbose",2);
    solver.initialize();

Solve the non-linear problem

The gsStaticNewton simply performs Newton-Raphson iterations under the hood. Solving the non-linear system is simply done by using:

    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsStatus status = solver.solve();
    GISMO_ASSERT(status==gsStatus::Success,"Solver failed");
    gsVector<> solVector = solver.solution();

Where the solution vector is obtained in the last line. As can be seen in the snipped above, the solver returns a status via gsStatus, which can be used to check whether the solver converged or not. Evaluation and export of the solutions is similar as in Tutorial: Linear Kirchhoff-Love shell analysis and Tutorial: Non-Linear Kirchhoff-Love shell analysis.

Annotated source file

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

 
#include <gismo.h>
 
#ifdef gsKLShell_ENABLED
#include <gsKLShell/gsKLShell.h>
#endif
 
#include <gsStructuralAnalysis/src/gsStaticSolvers/gsStaticNewton.h>
#include <gsStructuralAnalysis/src/gsStructuralAnalysisTools/gsStructuralAnalysisUtils.h>
 
using namespace gismo;
 
// Choose among various shell examples, default = Thin Plate
int main(int argc, char *argv[])
{
#ifdef gsKLShell_ENABLED
    bool plot  = false;
    index_t numRefine  = 1;
    index_t numElevate = 1;
 
    gsCmdLine cmd("Nonlinear static 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;
    gsReadFile<>("surfaces/paraboloid.xml", ori);
    ori.computeTopology();
 
    // p-refine
    if (numElevate!=0)
        ori.degreeElevate(numElevate);
    // h-refine
    for (int r =0; r < numRefine; ++r)
        ori.uniformRefine();
 
    gsMultiBasis<> bases(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
    bc.addCornerValue(boundary::southwest, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::southeast, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::northwest, 0.0, 0, -1); // (corner,value, patch, unknown)
    bc.addCornerValue(boundary::northeast, 0.0, 0, -1); // (corner,value, patch, unknown)
 
    gsPointLoads<real_t> pLoads = gsPointLoads<real_t>();
    gsVector<> point(2);
    gsVector<> load (3);
    point<< 0.5, 0.5 ; load << 0, 0.0, -1e4 ;
    pLoads.addLoad(point, load, 0 );
 
    // The surface force is defined in the physical space, i.e. 3D
    gsFunctionExpr<> force("0","0","0",3);
 
 
    // Define material parameters
    // The material parameters are defined in the physical domain as well (!)
    gsConstantFunction<> E (1.E9,3);
    gsConstantFunction<> nu(0.45,3);
    gsConstantFunction<> t (1E-2,3);
 
    // Define a linear material, see \ref gsMaterialMatrixLinear.h
    // The first parameter is the physical domain dimension
 
    // gsMaterialMatrixLinear<3,real_t> materialMatrix(ori,t,E,nu);
 
    std::vector<gsFunctionSet<>*> parameters{&E,&nu};
    gsOptionList options;
    options.addInt("Material","Material model: (0): SvK | (1): NH | (2): NH_ext | (3): MR | (4): Ogden",1);
    options.addInt("Implementation","Implementation: (0): Composites | (1): Analytical | (2): Generalized | (3): Spectral",1);
    gsMaterialMatrixBase<real_t>::uPtr materialMatrix = getMaterialMatrix<3,real_t>(ori,t,parameters,options);
 
    // Define the assembler.
    // The first parameter is the physical domain dimension
    // The second parameter is the real type
    // The third parameter controls the bending stiffness
    gsThinShellAssembler<3, real_t, true > assembler(ori,bases,bc,force,materialMatrix);
    assembler.setPointLoads(pLoads);
 
 
    // Function for the Jacobian
    gsStructuralAnalysisOps<real_t>::Jacobian_t Jacobian = [&assembler](gsVector<real_t> const &x, gsSparseMatrix<real_t> & m)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleMatrix(def);
        m = assembler.matrix();
        return status==ThinShellAssemblerStatus::Success;
    };
    // Function for the Residual
    gsStructuralAnalysisOps<real_t>::Residual_t Residual = [&assembler](gsVector<real_t> const &x, gsVector<real_t> & v)
    {
        gsMultiPatch<> def;
        assembler.constructSolution(x,def);
        ThinShellAssemblerStatus status = assembler.assembleVector(def);
        v = assembler.rhs();
        return status==ThinShellAssemblerStatus::Success;
    };
 
 
    assembler.assemble();
    gsSparseMatrix<> K = assembler.matrix();
    gsVector<> F = assembler.rhs();
 
    gsStaticNewton<real_t> solver(K,F,Jacobian,Residual);
    solver.options().setInt("verbose",2);
    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();
 
    gsMultiPatch<> displ = assembler.constructDisplacement(solVector);
    gsMultiPatch<> def = assembler.constructSolution(solVector);
 
    // Evaluate the solution on a reference point (parametric coordinates)
    // The reference points are stored column-wise
    gsMatrix<> refPoint(2,1);
    refPoint<<0.5,0.5;
    // 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", 1000, true);
 
        // Plot the membrane Von-Mises stresses on the geometry
        gsPiecewiseFunction<> VMm;
        assembler.constructStress(def,VMm,stress_type::von_mises_membrane);
        gsField<> membraneStress(def,VMm, true);
        gsWriteParaview(membraneStress,"MembraneStress",1000);
    }
    // ! [Export visualization in ParaView]
 
    return EXIT_SUCCESS;
#else
    GISMO_ERROR("The tutorial needs to be compiled with gsKLShell enabled");
    return EXIT_FAILED;
#endif
}// end main