gismo中用等几何解决线弹性问题的程序示例---未完待续2023.2.28

---

前言

gismo中用等几何解决线弹性问题

一、调用线弹性程序示例

1.1 对plateWithHole.xml文件的理解

算例来自文章：Isogeometric analysis: An overview and computer
implementation aspects

plateWithHole.xml文件

1.2 程序及注释

/// This is the 2D linear elasticity benchmark "Infinite plate with circular hole"
/// 这是二维线弹性基准“无限圆孔板”
/// as described in V.P.Nguyen, C.Anitescu, S.P.A.Bordas, T.Rabczuk, 2015
/// "Isogeometric analysis: An overview and computer implementation aspects".
/// 等几何分析:概述及计算机实现方面
/// Author: A.Shamanskiy (2016 - ...., TU Kaiserslautern)
#include <gismo.h>
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
#include<iostream>
using namespace std;

using namespace gismo;

int main(int argc, char* argv[]) {

    gsInfo << "This is the 2D linear elasticity benchmark: infinite plate with circular hole.\n";
    //gsInfo << " 这是二维线弹性问题:无限圆孔板 \n";

    //=====================================//
                // Input //
    //=====================================//

    std::string filename = "plateWithHole.xml";  // 打开    plateWithHole.xml文件
    index_t numUniRef = 1;     //  K细化次数 int型
    index_t numDegElev = 0;    //  升阶次数
    index_t numPlotPoints = 10000;    // preview软件画图的点数量
    bool plotMesh = false;  // preview画图时命令 136行：  是否画网格  false = 不画网格

    // 极简的终端用户界面
    gsCmdLine cmd("This is the 2D linear elasticity benchmark: infinite plate with circular hole.");
    //gsCmdLine cmd("这是二维线弹性问题:无限圆孔板");
    cmd.addInt("r", "refine", "Number of uniform refinement application", numUniRef);
    //cmd.addInt("r", "refine", " 细化次数    ", numUniRef);
    cmd.addInt("d", "degelev", "Number of degree elevation application：升阶次数？", numDegElev);
    // preview软件画图的点数量
    cmd.addInt("p", "points", "Number of points to plot to Paraview", numPlotPoints);
    cmd.addSwitch("m", "mesh", "Plot computational mesh：绘制计算网格", plotMesh);

    try { cmd.getValues(argc, argv); }   //  不太用看  不知道这个命令代表啥
    catch (int rv) { return rv; }   // 

    //=============================================//
        // Scanning geometry and creating bases:扫描几何和创建基函数 //
    //=============================================//

    // scanning geometry  扫描几何
    gsMultiPatch<> geometry;   //  定义一个多片
    gsReadFile<>(filename, geometry);    // 将plateWithHole.xml文件中的数据赋值给 geometry

    // 生成基函数
    gsMultiBasis<> basis(geometry);   // 生成基函数
    for (index_t i = 0; i < numDegElev; ++i)   // 升阶次数
        basis.degreeElevate();
    for (index_t i = 0; i < numUniRef; ++i)    // k细化（节点插入）次数
        basis.uniformRefine();
    
    //=============================================//
        // Setting loads and boundary conditions  设置载荷和边界条件 //
    //=============================================//
    // 应力的精确解
    gsFunctionExpr<> analyticalStresses("1-1/(x^2+y^2)*(3/2*cos(2*atan2(y,x)) + cos(4*atan2(y,x))) + 3/2/(x^2+y^2)^2*cos(4*atan2(y,x))",
        "-1/(x^2+y^2)*(1/2*cos(2*atan2(y,x)) - cos(4*atan2(y,x))) - 3/2/(x^2+y^2)^2*cos(4*atan2(y,x))",
        "-1/(x^2+y^2)*(1/2*sin(2*atan2(y,x)) + sin(4*atan2(y,x))) + 3/2/(x^2+y^2)^2*sin(4*atan2(y,x))", 2);
    // boundary load neumann BC 黎曼边界对应载荷边界条件  荷载BC   力的边界条件
    gsFunctionExpr<> traction("(-1+1/(x^2+y^2)*(3/2*cos(2*atan2(y,x)) + cos(4*atan2(y,x))) - 3/2/(x^2+y^2)^2*cos(4*atan2(y,x))) * (x==-4) +"
        "(-1/(x^2+y^2)*(1/2*sin(2*atan2(y,x)) + sin(4*atan2(y,x))) + 3/2/(x^2+y^2)^2*sin(4*atan2(y,x))) * (y==4)",
        "(1/(x^2+y^2)*(1/2*sin(2*atan2(y,x)) + sin(4*atan2(y,x))) - 3/2/(x^2+y^2)^2*sin(4*atan2(y,x))) * (x==-4) +"
        "(-1/(x^2+y^2)*(1/2*cos(2*atan2(y,x)) - cos(4*atan2(y,x))) - 3/2/(x^2+y^2)^2*cos(4*atan2(y,x))) * (y==4)", 2);
    // material parameters 材料参数
    real_t youngsModulus = 1.0e3;
    real_t poissonsRatio = 0.3;

    // boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件
    gsBoundaryConditions<> bcInfo;
    bcInfo.addCondition(0, boundary::north, condition_type::neumann, &traction);
    bcInfo.addCondition(0, boundary::west, condition_type::dirichlet, nullptr, 1); // last number is a component (coordinate) number
    bcInfo.addCondition(0, boundary::east, condition_type::dirichlet, nullptr, 0); // 最后一个数字是一个分量(坐标)的数字

    // source function, rhs  源函数，rhs
    gsConstantFunction<> g(0., 0., 2);   //  g = 0   0
    gsInfo << " g = \n " << g << "  \n ";

    //=============================================//
              // Assembling & solving  组装与解决 //
    //=============================================//

    // creating assembler  创建刚度矩阵
    gsElasticityAssembler<real_t> assembler(geometry, basis, bcInfo, g);
    assembler.options().setReal("YoungsModulus", youngsModulus);
    assembler.options().setReal("PoissonsRatio", poissonsRatio);
    gsInfo << "Assembling...\n";
    gsStopwatch clock;   //  计时
    clock.restart();    // 计时开始
    assembler.assemble();    //
   
    gsInfo << "Assembled a system (matrix and load vector) with "
        << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n";
    //  在clock.stop中装配assembler.numDofs自由度的系统(矩阵和负载矢量)
    gsInfo << "Solving...\n";
    clock.restart();  //  计时结束

#ifdef GISMO_WITH_PARDISO    //  
    gsSparseSolver<>::PardisoLLT solver(assembler.matrix());   // 将刚度矩阵传入solver中
    gsVector<> solVector = solver.solve(assembler.rhs());    // 将刚度矩阵传入solver中
    gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
    //用PardisoLDLT求解器解决了系统中的问题
#else
    gsSparseSolver<>::SimplicialLDLT solver(assembler.matrix());   //  assembler.matrix()传入solver
    gsVector<> solVector = solver.solve(assembler.rhs());   //  KU = p 的p传入solver
    gsInfo << "Solved the system with EigenLDLT solver in 在1中用EigenLDLT求解器求解系统 " << clock.stop() << "s.\n";
#endif

    //=============================================//
                      // Output //
    //=============================================//

    // 将位移构造为IGA函数
    gsMultiPatch<> solution;
    assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);
    // 构造应力张量
    gsPiecewiseFunction<> stresses;
    assembler.constructCauchyStresses(solution, stresses, stress_components::all_2D_vector);

    if (numPlotPoints > 0)
    {
        // 构造位移的IGA场(几何+解)
        gsField<> solutionField(assembler.patches(), solution);
        //  构造应力的IGA场(几何+解)
        gsField<> stressField(assembler.patches(), stresses, true);
        //  分析应力
        gsField<> analyticalStressField(assembler.patches(), analyticalStresses, false);
        // 创建一个容器，将所有字段绘制到一个Paraview文件中
        // 生成的Paraview文件：  文件名是什么   文件保存在哪
        std::map<std::string, const gsField<>*> fields;
        fields["Deformation"] = &solutionField;
        fields["Stress"] = &stressField;
        fields["StressAnalytical"] = &analyticalStressField;
        gsWriteParaviewMultiPhysics(fields, "plateWithHole", numPlotPoints, plotMesh);
        //gsInfo << "Open \"plateWithHole.pvd\" in Paraview for visualization. Stress wiggles on the left side are caused by "
          //  "a singularity in the parametrization.\n";
    }
    cout << "  \n ";
    gsInfo << "assembler.fixedDofs() = \n" << assembler.fixedDofs();
    cout << "  \n ";
    cout << "assembler.matrix() = \n" << assembler.matrix();
    cout << "  \n ";
    cout << "assembler.matrix().rows() = \n" << assembler.matrix().rows();
    cout << "  \n ";
    cout << "assembler.matrix().cols() = \n" << assembler.matrix().cols();
    cout << "  \n ";
    //gsInfo << "assembler.matrix() = \n" << assembler.matrix();

    return 0;
}

1.3 对边界力函数的理解

下边界应该是u = 0

总结 #pic_center

空格         空格

：

二维数

1

1

1