gismo程序示例:边长为 8 16 32 的长方体 受均布载荷


前言

只是为方便学习,不做其他用途,

一、

在这里插入图片描述

一、8*32面 受均布载荷

/// This is an example of using the linear elasticity solver on a 3D multi-patch geometry.
/// The problems is part of the EU project "Terrific".
///
///  示例:边长为 8   16   32  的长方体
///  8*32面  受均布载荷
/// Authors: O. Weeger (2012-1015, TU Kaiserslautern),
///          A.Shamanskiy (2016 - ...., TU Kaiserslautern)
#include <gismo.h>
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
#include <gsElasticity/gsGeoUtils.h>
#include<iostream>

using namespace std;
using namespace gismo;

int main(int argc, char* argv[])
{
    gsInfo << "Testing the linear elasticity solver in 3D-线弹性求解器的三维测试.\n";

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

    //std::string filename("terrific.xml");//初始数据文件
    std::string filename("test.xml");//初始数据文件
    real_t youngsModulus = 1e5;//杨氏模量
    real_t poissonsRatio = 0.3;//泊松比
    index_t numUniRef = 0;//节点插入数
    index_t numDegElev = 1;//升阶次数
    index_t numPlotPoints = 10000;//preview软件画图的点数量


    // minimalistic user interface for terminal 终端最简用户界面
    gsCmdLine cmd("Testing the linear elasticity solver in 3D.");// 定义一个gsCmdLine类 命名为cmd
    cmd.addInt("r", "refine", "Number of uniform refinement application", numUniRef);
    cmd.addInt("d", "degelev", "Number of degree elevation application", numDegElev);
    cmd.addInt("p", "points", "Number of points to plot to Paraview", numPlotPoints);
    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
    // creating basis 生成基函数
    gsMultiBasis<> basis(geometry);
    for (index_t i = 0; i < numDegElev; ++i) // 升阶次数
        basis.degreeElevate();
    for (index_t i = 0; i < numUniRef; ++i) // k细化(节点插入)次数
        basis.uniformRefine();
    gsInfo << basis;

    //=====================================================================//
        // Setting loads and boundary conditions  设置载荷和边界条件 //
    //=====================================================================//

    // source function, rhs  源函数?-解析解?
    gsConstantFunction<> f(0., 0., 0., 3);   //  g = 0   0   0
    gsInfo << " f = \n " << f << "  \n ";
    // surface load, neumann BC 黎曼边界对应载荷边界条件  荷载BC   力的边界条件
    int W = 8;
    int H = 32;
    int g1 = 1e3 / (W * H);
    gsConstantFunction<> g(0, 0, g1, 3);   //  g = 0   0
    //gsConstantFunction<> g(0, 0, 100, 3);   //  g = 0   0
    gsInfo << " g = \n " << g << "  \n ";

    // boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件
    gsBoundaryConditions<> bcInfo;
    // Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BC
    for (index_t d = 0; d < 3; d++)
    {
        bcInfo.addCondition(0, boundary::west, condition_type::dirichlet, 0, d);
    }
    // Neumann BC are imposed as one function 将 Neumann BC 作为一个函数
    bcInfo.addCondition(0, boundary::east, condition_type::neumann, &g);


    //=====================================================================//
              // Assembling & solving //
    //=====================================================================//

    // creating assembler   创建刚度矩阵?
    gsElasticityAssembler<real_t> assembler(geometry, basis, bcInfo, f);
    assembler.options().setReal("YoungsModulus", youngsModulus);
    assembler.options().setReal("PoissonsRatio", poissonsRatio);
    //assembler.options().setInt("DirichletValues", dirichlet::l2Projection);
    gsInfo << "Assembling...\n";
    gsStopwatch clock;
    clock.restart();
    assembler.assemble();
    //assembler.matrix();//总刚
    gsInfo << "Assembled a system with "
        << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n";

    gsInfo << "Solving...\n";
    clock.restart();

#ifdef GISMO_WITH_PARDISO
    gsSparseSolver<>::PardisoLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#else
    gsSparseSolver<>::SimplicialLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with EigenLDLT solver in " << clock.stop() << "s.\n";
#endif

    // 输出的位移和总刚都是将   位移不变的自由度去掉得到的结果
    /*gsInfo << " 位移:solVector = \n" << solVector;
    cout << " \n ";*/
    double max = solVector.maxCoeff();
    cout << "\n 最大节点位移 solVector = \n" << max << endl;
    double Solution = 6.667e-5;
    double error = (max - Solution) / max;
    cout << "\n 误差 = " << abs(error) * 100 << "%" << endl;

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

    // constructing solution as an IGA function
    gsMultiPatch<> solution;
    assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);
    /*gsInfo << " solution = \n" << solution;
    cout << " \n ";*/
    // constructing stresses
    gsPiecewiseFunction<> stresses;
    assembler.constructCauchyStresses(solution, stresses, stress_components::von_mises);

    if (numPlotPoints > 0)
    {
        // constructing an IGA field (geometry + solution)
        gsField<> solutionField(assembler.patches(), solution);
        gsField<> stressField(assembler.patches(), stresses, true);
        // creating a container to plot all fields to one Paraview file
        std::map<std::string, const gsField<>*> fields;
        fields["Deformation"] = &solutionField;
        fields["von Mises"] = &stressField;
        gsWriteParaviewMultiPhysics(fields, "test3_le", numPlotPoints);
        gsInfo << "Open \"test3_le.pvd\" in Paraview for visualization.\n";
    }

    return 0;
}

二、最小的面(8*16)受均布载荷

/// This is an example of using the linear elasticity solver on a 3D multi-patch geometry.
/// The problems is part of the EU project "Terrific".
///
///  示例:边长为 8   16   32  的长方体
///  最小的面(8*16)受均布载荷
/// Authors: O. Weeger (2012-1015, TU Kaiserslautern),
///          A.Shamanskiy (2016 - ...., TU Kaiserslautern)
#include <gismo.h>
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
#include <gsElasticity/gsGeoUtils.h>
#include<iostream>


using namespace std;
using namespace gismo;

int main(int argc, char* argv[])
{
    gsInfo << "Testing the linear elasticity solver in 3D-线弹性求解器的三维测试.\n";

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

    //std::string filename("terrific.xml");//初始数据文件
    std::string filename("test.xml");//初始数据文件
    real_t youngsModulus = 74e9;//杨氏模量
    real_t poissonsRatio = 0.33;//泊松比
    index_t numUniRef = 2;//节点插入数
    index_t numDegElev = 1;//升阶次数
    index_t numPlotPoints = 10000;//preview软件画图的点数量


    // minimalistic user interface for terminal 终端最简用户界面
    gsCmdLine cmd("Testing the linear elasticity solver in 3D.");// 定义一个gsCmdLine类 命名为cmd
    cmd.addInt("r", "refine", "Number of uniform refinement application", numUniRef);
    cmd.addInt("d", "degelev", "Number of degree elevation application", numDegElev);
    cmd.addInt("p", "points", "Number of points to plot to Paraview", numPlotPoints);
    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
    // creating basis 生成基函数
    gsMultiBasis<> basis(geometry);
    for (index_t i = 0; i < numDegElev; ++i) // 升阶次数
        basis.degreeElevate();
    for (index_t i = 0; i < numUniRef; ++i) // k细化(节点插入)次数
        basis.uniformRefine();
    gsInfo << basis;

    //=====================================================================//
        // Setting loads and boundary conditions  设置载荷和边界条件 //
    //=====================================================================//

    // source function, rhs  源函数?-解析解?
    gsConstantFunction<> f(0., 0., 0., 3);   //  g = 0   0   0
    gsInfo << " f = \n " << f << "  \n ";
    // surface load, neumann BC 黎曼边界对应载荷边界条件  荷载BC   力的边界条件
    int W = 8;
    int H = 16;
    int g1 = 2e7 / (W * H);
    gsConstantFunction<> g(g1, 0, 0, 3);   //  g = 0   0
    gsInfo << " g = \n " << g << "  \n ";

    // boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件
    gsBoundaryConditions<> bcInfo;
    // Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BC
    for (index_t d = 0; d < 3; d++)
    {
        bcInfo.addCondition(0, boundary::south, condition_type::dirichlet, 0, d);
    }
    // Neumann BC are imposed as one function 将 Neumann BC 作为一个函数
    bcInfo.addCondition(0, boundary::north, condition_type::neumann, &g);


    //=====================================================================//
              // Assembling & solving //
    //=====================================================================//

    // creating assembler   创建刚度矩阵?
    gsElasticityAssembler<real_t> assembler(geometry, basis, bcInfo, f);
    assembler.options().setReal("YoungsModulus", youngsModulus);
    assembler.options().setReal("PoissonsRatio", poissonsRatio);
    //assembler.options().setInt("DirichletValues", dirichlet::l2Projection);
    gsInfo << "Assembling...\n";
    gsStopwatch clock;
    clock.restart();
    assembler.assemble();
    gsInfo << "Assembled a system with "
        << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n";

    gsInfo << "Solving...\n";
    clock.restart();

#ifdef GISMO_WITH_PARDISO
    gsSparseSolver<>::PardisoLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#else
    gsSparseSolver<>::SimplicialLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with EigenLDLT solver in " << clock.stop() << "s.\n";
#endif

    // 输出的位移和总刚都是将   位移不变的自由度去掉得到的结果
    /*gsInfo << " 位移:solVector = \n" << solVector;
    cout << " \n ";*/
    double max = solVector.maxCoeff();
    cout << "\n 最大节点位移 solVector = \n" << max << endl;
    double Solution = 6.667e-5;
    double error = (max - Solution) / max;
    cout << "\n 误差 = " << abs(error) * 100 << "%" << endl;

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

    // constructing solution as an IGA function
    gsMultiPatch<> solution;
    assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);
    /*gsInfo << " solution = \n" << solution;
    cout << " \n ";*/
    // constructing stresses
    gsPiecewiseFunction<> stresses;
    assembler.constructCauchyStresses(solution, stresses, stress_components::von_mises);

    if (numPlotPoints > 0)
    {
        // constructing an IGA field (geometry + solution)
        gsField<> solutionField(assembler.patches(), solution);
        gsField<> stressField(assembler.patches(), stresses, true);
        // creating a container to plot all fields to one Paraview file
        std::map<std::string, const gsField<>*> fields;
        fields["Deformation"] = &solutionField;
        fields["von Mises"] = &stressField;
        gsWriteParaviewMultiPhysics(fields, "test2_le", numPlotPoints);
        gsInfo << "Open \"test2_le.pvd\" in Paraview for visualization.\n";
    }

    return 0;
}

三、最大的面受均布载荷

/// This is an example of using the linear elasticity solver on a 3D multi-patch geometry.
/// The problems is part of the EU project "Terrific".
///
/// 最大的面受均布载荷
/// Authors: O. Weeger (2012-1015, TU Kaiserslautern),
///          A.Shamanskiy (2016 - ...., TU Kaiserslautern)
#include <gismo.h>
#include <gsElasticity/gsElasticityAssembler.h>
#include <gsElasticity/gsWriteParaviewMultiPhysics.h>
#include <gsElasticity/gsGeoUtils.h>
#include<iostream>

using namespace std;
using namespace gismo;

int main(int argc, char* argv[])
{
    gsInfo << "Testing the linear elasticity solver in 3D-线弹性求解器的三维测试.\n";

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

    //std::string filename("terrific.xml");//初始数据文件
    std::string filename("test.xml");//初始数据文件
    real_t youngsModulus = 74e9;//杨氏模量
    real_t poissonsRatio = 0.33;//泊松比
    index_t numUniRef = 2;//节点插入数
    index_t numDegElev = 1;//升阶次数
    index_t numPlotPoints = 10000;//preview软件画图的点数量


    // minimalistic user interface for terminal 终端最简用户界面
    gsCmdLine cmd("Testing the linear elasticity solver in 3D.");// 定义一个gsCmdLine类 命名为cmd
    cmd.addInt("r", "refine", "Number of uniform refinement application", numUniRef);
    cmd.addInt("d", "degelev", "Number of degree elevation application", numDegElev);
    cmd.addInt("p", "points", "Number of points to plot to Paraview", numPlotPoints);
    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
    // creating basis 生成基函数
    gsMultiBasis<> basis(geometry);
    for (index_t i = 0; i < numDegElev; ++i) // 升阶次数
        basis.degreeElevate();
    for (index_t i = 0; i < numUniRef; ++i) // k细化(节点插入)次数
        basis.uniformRefine();
    gsInfo << basis;

    //=====================================================================//
        // Setting loads and boundary conditions  设置载荷和边界条件 //
    //=====================================================================//

    // source function, rhs  源函数?-解析解?
    gsConstantFunction<> f(0., 0., 0., 3);   //  g = 0   0   0
    gsInfo << " f = \n " << f << "  \n ";
    // surface load, neumann BC 黎曼边界对应载荷边界条件  荷载BC   力的边界条件
    int W = 16;
    int H = 32;
    int g1 = 2e7 / (W * H);
    gsConstantFunction<> g(0, g1, 0, 3); //  g = 0   0
    gsInfo << " g = \n " << g << "  \n ";

    // boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件
    gsBoundaryConditions<> bcInfo;
    // Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BC
    for (index_t d = 0; d < 3; d++)
    {
        bcInfo.addCondition(0, boundary::back, condition_type::dirichlet, 0, d);
    }
    // Neumann BC are imposed as one function 将 Neumann BC 作为一个函数
    bcInfo.addCondition(0, boundary::front, condition_type::neumann, &g);


    //=====================================================================//
              // Assembling & solving //
    //=====================================================================//

    // creating assembler   创建刚度矩阵?
    gsElasticityAssembler<real_t> assembler(geometry, basis, bcInfo, f);
    assembler.options().setReal("YoungsModulus", youngsModulus);
    assembler.options().setReal("PoissonsRatio", poissonsRatio);
    assembler.options().setInt("DirichletValues", dirichlet::l2Projection);
    gsInfo << "Assembling...\n";
    gsStopwatch clock;
    clock.restart();
    assembler.assemble();
    gsInfo << "Assembled a system with "
        << assembler.numDofs() << " dofs in " << clock.stop() << "s.\n";

    gsInfo << "Solving...\n";
    clock.restart();

#ifdef GISMO_WITH_PARDISO
    gsSparseSolver<>::PardisoLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#else
    gsSparseSolver<>::SimplicialLDLT solver(assembler.matrix());
    gsVector<> solVector = solver.solve(assembler.rhs());
    gsInfo << "Solved the system with EigenLDLT solver in " << clock.stop() << "s.\n";
#endif

    // 输出的位移和总刚都是将   位移不变的自由度去掉得到的结果
    //gsInfo << " 位移:solVector = \n" << solVector;
    //cout << " \n ";
    double max = solVector.maxCoeff();
    cout << "\n 最大节点位移 solVector = \n" << max << endl;
    double Solution = 4.696e-6;
    double error = (max - Solution) / max;
    cout << "\n 误差 = " << abs(error) * 100 << "%" << endl;

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

    // constructing solution as an IGA function
    gsMultiPatch<> solution;
    assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);
    gsInfo << " solution = \n" << solution;
    cout << " \n ";
    // constructing stresses
    gsPiecewiseFunction<> stresses;
    assembler.constructCauchyStresses(solution, stresses, stress_components::von_mises);

    if (numPlotPoints > 0)
    {
        // constructing an IGA field (geometry + solution)
        gsField<> solutionField(assembler.patches(), solution);
        gsField<> stressField(assembler.patches(), stresses, true);
        // creating a container to plot all fields to one Paraview file
        std::map<std::string, const gsField<>*> fields;
        fields["Deformation"] = &solutionField;
        fields["von Mises"] = &stressField;
        gsWriteParaviewMultiPhysics(fields, "test1_le", numPlotPoints);
        gsInfo << "Open \"test1_le.pvd\" in Paraview for visualization.\n";
    }

    return 0;
}
  • 0
    点赞
  • 0
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值