文章目录
- 前言
- 一、
- 一、8*32面 受均布载荷
- 二、最小的面（8*16）受均布载荷
- 三、最大的面受均布载荷
前言
只是为方便学习，不做其他用途，
一、
一、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;
}
