前言

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

一、

在这里插入图片描述

一、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类 命名为cmdcmd.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   0gsInfo << " 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   0gsInfo << " g = \n " << g << "  \n ";// boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件gsBoundaryConditions<> bcInfo;// Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BCfor (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_PARDISOgsSparseSolver<>::PardisoLDLT solver(assembler.matrix());gsVector<> solVector = solver.solve(assembler.rhs());gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#elsegsSparseSolver<>::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 functiongsMultiPatch<> solution;assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);/*gsInfo << " solution = \n" << solution;cout << " \n ";*/// constructing stressesgsPiecewiseFunction<> 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 filestd::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类 命名为cmdcmd.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   0gsInfo << " 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   0gsInfo << " g = \n " << g << "  \n ";// boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件gsBoundaryConditions<> bcInfo;// Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BCfor (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_PARDISOgsSparseSolver<>::PardisoLDLT solver(assembler.matrix());gsVector<> solVector = solver.solve(assembler.rhs());gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#elsegsSparseSolver<>::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 functiongsMultiPatch<> solution;assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);/*gsInfo << " solution = \n" << solution;cout << " \n ";*/// constructing stressesgsPiecewiseFunction<> 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 filestd::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类 命名为cmdcmd.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   0gsInfo << " 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   0gsInfo << " g = \n " << g << "  \n ";// boundary conditions 边界条件  黎曼边界对应载荷边界条件    dirichlete对应位移边界条件gsBoundaryConditions<> bcInfo;// Dirichlet BC are imposed separately for every component (coordinate) 对每个分量(坐标)分别施加 Dirichlet BCfor (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_PARDISOgsSparseSolver<>::PardisoLDLT solver(assembler.matrix());gsVector<> solVector = solver.solve(assembler.rhs());gsInfo << "Solved the system with PardisoLDLT solver in " << clock.stop() << "s.\n";
#elsegsSparseSolver<>::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 functiongsMultiPatch<> solution;assembler.constructSolution(solVector, assembler.allFixedDofs(), solution);gsInfo << " solution = \n" << solution;cout << " \n ";// constructing stressesgsPiecewiseFunction<> 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 filestd::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;
}