Tutorial: Linear Kirchhoff-Love shell analysis

This example solves a linear shell problem using the gsThinShellAssembler. The goals of the tutorial are the following:

• Construct a gsMaterialMatrixBase object, given material parameters
• Construct a gsThinShellAssembler object, given boundary conditions
• Solve a linear problem and export the solutions

Preliminaries

Like every example in G+Smo, we start our file with some includes. In this case, we include the gsThinShellAssembler.h file containing the assembler for the Kirchhoff-Love shell, and the getMaterialMatrix.h used for getting the material matrix.

#include <gsKLShell/gsKLShell.h>

The command line options for this file are relatively easy. One can control the number of mesh refinements and degree elevations, and specify a plot flag if needed.

    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; }

Geometry and basis

The problem at hand will be the Scordelis-Lo Roof problem from the shell obstactle course, see this paper. The geometry is stored in a file, hence it is read into the original geometry ori via gsReadFile.

    // Initialize [ori]ginal and [def]ormed geometry
    gsMultiPatch<> ori;
    std::string fileName = "surfaces/scordelis_lo_roof.xml";
    gsReadFile<>(fileName, ori);

For the analysis, we will use the geometry defined in the basis that we want for analysis. This is not essential, but it makes it easier to construct the geometry later using the gsThinShellAssembler class

    // p-refine
    if (numElevate!=0)
        ori.degreeElevate(numElevate);
    // h-refine
    for (int r =0; r < numRefine; ++r)
        ori.uniformRefine();

The basis is simply obtained using gsMultiBasis

    gsMultiBasis<> bases(ori);

Boundary conditions and loads

The boundary conditions of the problem are assigned via gsBoundaryConditions. On the east and west boundary of the geometry, the \(y\) and \(z\) displacements are fixed, and to make the problem well-defined, the south-west corner is fixed in \(x\). The boundary conditions are defined by first specifying the side (boxSide), then the type of condition (condition_type), then the function (nullptr here), the unknown (always 0, since displacements are the only unknown field), whether the condition is applied in the parametric domain (always false) and finally the component of the solution field (-1 for all, 0 for \(x\), 1 for \(y\) and 2 for \(z\)):

    // 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.addCondition(boundary::west, condition_type::dirichlet, 0, 0 ,false, 1 );
    bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0 ,false, 2 );
 
    bc.addCornerValue(boundary::southwest, 0.0, 0, 0); // (corner,value, patch, unknown)
 
    bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false, 1 );
    bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false, 2 );

The surface load on the shell is provided by a gsFunctionSet, which in this case is a derived gsFunctionExpr.

    // The surface force is defined in the physical space, i.e. 3D
    gsFunctionExpr<> force("0","0","-90",3);

Material definition

The material law used in gsThinShellAssembler is specified using a separate class, derived from gsMaterialMatrixBase. The family of this class includes several material models, and can be evaluated in every point (including a thickness coordinate) using gsMaterialMatrixEval or it can be integrated using gsMaterialMatrixIntegrate. Depending on the options specified to the latter classes, different constitutive quantities can be obtaied, including the normal force and bending moment.

To define a meterial matrix, we can use the getMaterialMatrix function. This functions needs the undeformed geometry, the thickness, a set of parameters and some options:

    // Define material parameters
    // The material parameters are defined in the physical domain as well (!)
    gsConstantFunction<> E (4.32E8,3);
    gsConstantFunction<> nu(0.0   ,3);
    gsConstantFunction<> t (0.25  ,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",0);
    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);

In case of "Material" 0, a gsMaterialMatrixLinear is provided, which only takes the Young's Modulus ( \(E\)) and the Poisson's ratio ( \(\nu\)). For the linear material matrix, thickness integration is simply done by multiplying with the thickness moments. For hyperelastic material models (see gsMaterialMatrixNonlinear), more parameters can be required, and through-thickness integration is done numerically. We refer to the page The gsMaterialMatrixBase for more details on material matrices.

Assembling and solving the linear problem

When the geometry, the basis, the boundary conditions and forces and the material are specified, the gsThinShellAssembler is simply constructed as follows:

    // 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);

Here, the first template parameter is the physical dimension, the second is the real type and the third is a boolean specifying whether bending stiffness contributions should be assembled.

Assemlbly of the linear stiffness matrix and external load vector is simply performed by calling

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

Having the stiffness matrix and force vector assembled, the linear shell problem can be solved accordingly. This is done with a simple gsSparseSolver:

    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsVector<> solVector;
    gsSparseSolver<>::CGDiagonal solver;
    solver.compute( matrix );
    solVector = solver.solve(vector);

Exporting the solution

Solving the linear problem gives the displacements of the control points of the geometry. To construct the displacement field, these control points form a solution field using the same basis as used for analysis. To construct the deformed geometry, their values need to be added to the control points of the original geometry. For ease of use, these functions are available in the gsThinShellAssembler:

    gsMultiPatch<> displ = assembler.constructDisplacement(solVector);
    gsMultiPatch<> def = assembler.constructSolution(solVector);

Using the displacement field and the deformed geometry, one can make some evaluations like with any other gsFunction, gsGeometry or gsMultiPatch, for example to print the displacements in a point:

    // Evaluate the solution on a reference point (parametric coordinates)
    // The reference points are stored column-wise
    gsMatrix<> refPoint(2,1);
    refPoint<<0.5,1;
    // 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";

In addition, the solutions can be exported to Paraview using gsField:

    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);
    }

Annotated source file

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

 
#include <gismo.h>
 
#include <gsKLShell/gsKLShell.h>
 
using namespace gismo;
 
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 = "surfaces/scordelis_lo_roof.xml";
    gsReadFile<>(fileName, ori);
 
    // 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.addCondition(boundary::west, condition_type::dirichlet, 0, 0 ,false, 1 );
    bc.addCondition(boundary::west, condition_type::dirichlet, 0, 0 ,false, 2 );
 
    bc.addCornerValue(boundary::southwest, 0.0, 0, 0); // (corner,value, patch, unknown)
 
    bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false, 1 );
    bc.addCondition(boundary::east, condition_type::dirichlet, 0, 0 ,false, 2 );
 
    // The surface force is defined in the physical space, i.e. 3D
    gsFunctionExpr<> force("0","0","-90",3);
 
 
    // Define material parameters
    // The material parameters are defined in the physical domain as well (!)
    gsConstantFunction<> E (4.32E8,3);
    gsConstantFunction<> nu(0.0   ,3);
    gsConstantFunction<> t (0.25  ,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",0);
    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.assemble();
    gsSparseMatrix<> matrix = assembler.matrix();
    gsVector<> vector = assembler.rhs();
 
    gsInfo<<"Solving system with "<<assembler.numDofs()<<" DoFs\n";
    gsVector<> solVector;
    gsSparseSolver<>::CGDiagonal solver;
    solver.compute( matrix );
    solVector = solver.solve(vector);
 
    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,1;
    // 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;
 
}// end main