Tutorial: Linear Kirchhoff-Love shell analysis on multi-patches using gsUnstructuredSplines

This example solves a linear shell problem on a multi-patch using the gsThinShellAssembler and the Unstructured splines module module. The goals of the tutorial are the following:

• Define an unstructured spline
• Solve the linear shell equation on the multi-patch basis

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

The problem to be solved corresponds to the hyperboloid shell example from [1].

Geometry and basis

Firstly, the geometry is loaded in the same way as we load a geometry in Tutorial: Linear Kirchhoff-Love shell analysis. To be sure that our geometry has interfaces and boundaries, we compute its topology:

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

Refinement and elevation of the geometry (and the basis) are performed on the multi-patch geometry as is. This is the same as in Tutorial: Linear Kirchhoff-Love shell analysis.

To perform isogeometric analysis on the multi-patch geometry and a smooth corresponding multi-basis, we will use the gsMappedBasis and gsMappedSpline. These objects use a locally defined (non-smooth) basis and a sparse mapping matrix. The smooth basis is defined using the D-Patch construction from [2], available in G+Smo using gsDPatch. The D-Patch is constructed on the gsMultiBasis belonging to our original geometry.

The construction of the smooth basis and the corresponding geometry is done as follows:

    // Construct a multi-basis object
    gsMultiBasis<> bases(ori);
    // Compute the smooth D-Patch construction
    gsDPatch<2,real_t> dpatch(bases);
    dpatch.options().setSwitch("SharpCorners",false);
    dpatch.compute();
 
    // Store the basis mapping matrix
    gsSparseMatrix<> global2local;
    dpatch.matrix_into(global2local);
    global2local = global2local.transpose();
    // Overwrite the local basis (might change due to D-Patch)
    bases = dpatch.localBasis();
 
    // Store the modified, smooth geometry
    tmp = ori; // temporary storage
    ori = dpatch.exportToPatches(tmp);
 
    // Construct a mapped basis, using the local basis and the mapper
    gsMappedBasis<2,real_t> mbasis;
    mbasis.init(bases,global2local);

Boundary conditions and loads

The boundary conditions and external forces correspond to the hyperboloid shell case from [1], i.e. one edge of the hyperboloid is fully fixed. The boundary conditions are consequently applied to patch 0 and 1. See Boundary conditions and loads for more details on applying boundary conditions.

    // 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(0,boundary::west, condition_type::dirichlet, 0, 0 ,false, -1 );
    bc.addCondition(0,boundary::west, condition_type::clamped  , 0, 0 ,false, -1 );
    bc.addCondition(1,boundary::west, condition_type::dirichlet, 0, 0 ,false, -1 );
    bc.addCondition(1,boundary::west, condition_type::clamped  , 0, 0 ,false, -1 );

Similarly, we define a surface force corresponding to the benchmark problem

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

Configuration of the assembler

The configuration of the material matrix and the shell assembler is the same as in Tutorial: Linear Kirchhoff-Love shell analysis, just with different material parameters, corresponding to the benchmark problem in [1].

Solving the problem

Assembly of the problem, as well as solving it, is the same as in Tutorial: Linear Kirchhoff-Love shell analysis.

Exporting the solution

Construction of the solution works slightly different when a gsMappedBasis is used for analysis. Instead of storing the solutions in a gsMultiPatch, we now store them in a gsMappedSpline. To define the deformed geometry on top of the original geometry, one can use the gsFunctionSum class, which basically evaluates two functions on a parametric point, and sums their values. The full code for the construction of the deformed geometry is here:

    // We will construct a mapped spline using all coefficients
    gsMatrix<real_t> solFull = assembler.fullSolutionVector(solVector);
 
    // 2. Reshape all the coefficients to a Nx3 matrix
    GISMO_ASSERT(solFull.rows() % 3==0,"Rows of the solution vector does not match the number of control points");
    solFull.resize(solFull.rows()/3,3);
 
    // 3. Make the mapped spline
    gsMappedSpline<2,real_t> mspline(mbasis,solFull);
 
    // Deformation is, on every parametric point, the original geometry + the mapped spline
    gsFunctionSum<real_t> def(&ori,&mspline);

As can be seen in the following snippet, the constructed solution object can be used for evaluation:

    // Evaluate the solution on a reference point (parametric coordinates)
    // The reference points are stored column-wise
    gsMatrix<> refPoint(2,1);
    index_t refPatch;
    // In this case, the reference point is on the east boundary of patch 4
    refPoint<<1.0,0.5;
    refPatch = 4;
    // Compute the values
    gsVector<> physpoint = def  .piece(refPatch).eval(refPoint);
    gsVector<> refVals = mspline.piece(refPatch).eval(refPoint);
    // gsInfo << "Displacement at reference point: "<<numVal<<"\n";
    gsInfo  << "Displacement at reference point ("
            <<physpoint.transpose()<<"): "
            <<": "<<refVals.transpose()<<"\n";
    gsInfo  <<"The elastic energy is: "<<0.5*solVector.transpose()*matrix*solVector<<"\n";

And for visualization:

    if (plot)
    {
        // Plot the displacements on the deformed geometry
        gsField<> solField(ori, mspline,true);
        gsInfo<<"Plotting in Paraview...\n";
        gsWriteParaview<>( solField, "Deformation", 1000, true);
 
        // Plot the membrane Von-Mises stresses on the geometry
        gsPiecewiseFunction<> VMm;
        assembler.constructStress(ori,def,VMm,stress_type::von_mises_membrane);
        gsField<> membraneStress(ori,VMm, true);
        gsWriteParaview(membraneStress,"MembraneStress",1000);
    }

References

[1] Bathe, Klaus-Jürgen, Alexander Iosilevich, and Dominique Chapelle. "An evaluation of the MITC shell elements." Computers & Structures 75.1 (2000): 1-30.

[2] Toshniwal D, Speleers H, Hughes TJ. Smooth cubic spline spaces on unstructured quadrilateral meshes with particular emphasis on extraordinary points: Geometric design and isogeometric analysis considerations. Computer Methods in Applied Mechanics and Engineering. 2017 Dec 1;327:411-58.

Annotated source file

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

 
#include <gismo.h>
 
#include <gsUnstructuredSplines/src/gsDPatch.h>
 
#ifdef gsKLShell_ENABLED
#include <gsKLShell/src/gsThinShellAssembler.h>
#include <gsKLShell/src/getMaterialMatrix.h>
#include <gsKLShell/src/gsFunctionSum.h>
#endif
 
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("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, tmp;
    std::string fileName = "surfaces/6p_hyperboloid.xml";
    gsReadFile<>(fileName, ori);
    ori.computeTopology();
 
    // p-refine
    if (numElevate!=0)
        ori.degreeElevate(numElevate);
    // h-refine
    for (int r =0; r < numRefine; ++r)
        ori.uniformRefine();
 
    // Construct a multi-basis object
    gsMultiBasis<> bases(ori);
    // Compute the smooth D-Patch construction
    gsDPatch<2,real_t> dpatch(bases);
    dpatch.options().setSwitch("SharpCorners",false);
    dpatch.compute();
 
    // Store the basis mapping matrix
    gsSparseMatrix<> global2local;
    dpatch.matrix_into(global2local);
    global2local = global2local.transpose();
    // Overwrite the local basis (might change due to D-Patch)
    bases = dpatch.localBasis();
 
    // Store the modified, smooth geometry
    tmp = ori; // temporary storage
    ori = dpatch.exportToPatches(tmp);
 
    // Construct a mapped basis, using the local basis and the mapper
    gsMappedBasis<2,real_t> mbasis;
    mbasis.init(bases,global2local);
 
    // 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(0,boundary::west, condition_type::dirichlet, 0, 0 ,false, -1 );
    bc.addCondition(0,boundary::west, condition_type::clamped  , 0, 0 ,false, -1 );
    bc.addCondition(1,boundary::west, condition_type::dirichlet, 0, 0 ,false, -1 );
    bc.addCondition(1,boundary::west, condition_type::clamped  , 0, 0 ,false, -1 );
 
    // The surface force is defined in the physical space, i.e. 3D
    // The load is -8000*thickness
    gsFunctionExpr<> force("0","0","-80",3);
 
    // Define material parameters
    // The material parameters are defined in the physical domain as well (!)
    gsConstantFunction<> E (2E11,3);
    gsConstantFunction<> nu(0.3 ,3);
    gsConstantFunction<> t (0.01,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.setSpaceBasis(mbasis);
 
    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);
 
    // We will construct a mapped spline using all coefficients
    gsMatrix<real_t> solFull = assembler.fullSolutionVector(solVector);
 
    // 2. Reshape all the coefficients to a Nx3 matrix
    GISMO_ASSERT(solFull.rows() % 3==0,"Rows of the solution vector does not match the number of control points");
    solFull.resize(solFull.rows()/3,3);
 
    // 3. Make the mapped spline
    gsMappedSpline<2,real_t> mspline(mbasis,solFull);
 
    // Deformation is, on every parametric point, the original geometry + the mapped spline
    gsFunctionSum<real_t> def(&ori,&mspline);
 
 
    // Evaluate the solution on a reference point (parametric coordinates)
    // The reference points are stored column-wise
    gsMatrix<> refPoint(2,1);
    index_t refPatch;
    // In this case, the reference point is on the east boundary of patch 4
    refPoint<<1.0,0.5;
    refPatch = 4;
    // Compute the values
    gsVector<> physpoint = def  .piece(refPatch).eval(refPoint);
    gsVector<> refVals = mspline.piece(refPatch).eval(refPoint);
    // gsInfo << "Displacement at reference point: "<<numVal<<"\n";
    gsInfo  << "Displacement at reference point ("
            <<physpoint.transpose()<<"): "
            <<": "<<refVals.transpose()<<"\n";
    gsInfo  <<"The elastic energy is: "<<0.5*solVector.transpose()*matrix*solVector<<"\n";
 
    // ! [Export visualization in ParaView]
    if (plot)
    {
        // Plot the displacements on the deformed geometry
        gsField<> solField(ori, mspline,true);
        gsInfo<<"Plotting in Paraview...\n";
        gsWriteParaview<>( solField, "Deformation", 1000, true);
 
        // Plot the membrane Von-Mises stresses on the geometry
        gsPiecewiseFunction<> VMm;
        assembler.constructStress(ori,def,VMm,stress_type::von_mises_membrane);
        gsField<> membraneStress(ori,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_FAILURE;
#endif
}// end main