Intro 1
第一个例子演示了怎样创建一个Mesh Object。
头文件：1
2
3
using namespace libMesh; //所有的东西都在libMesh命名空间中
入口函数代码：1
LibMeshInit init (argc, argv); // 初始化，这是必要的一步，因为libMesh可能依赖其他的比如MPI、PETSc等函数库
因为该程序的用法是：1
./ex1 -d DIM input_mesh_name [-o output_mesh_name]
所有这里会有一个参数个数的判断，如果少于4个参数，就会报错，提示正确的用法：1
libmesh_error_msg("Usage: " << argv[0] << " -d 2 in.mesh [-o out.mesh]");
以下是对mesh object的操作：1
2
3
4Mesh mesh(init.comm()); // 在默认通信子上创建一个Mesh
mesh.read(argv[3]); //读入网格文件
mesh.print_info(); //在屏幕上输出网格信息
mesh.write(argv[5]); //输出网格文件
Intro 2
第二个例子演示了怎样创建一个标量方程。同时引入了PETSc，默认情形下方程数据都存入PETSc vector中。
头文件：1
2
3
4
5
6
7
8
9
10
11
12
//实现网格功能的基本头文件
// 定义简单的稳态系统
using namespace libMesh;
源码：1
2Mesh mesh(init.comm()); // 创建mesh
MeshTools::Genertion::build_square (mesh,5,5); // 划分成5×5的方形grid
这样创建的网格就有55=25个单元elements，66=36个节点nodes。1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17EquationSystems eqution_systems (mesh); // 创建方程，需要将mesh传递给equation
equation_systems.add_system<TransientLinearImplicitSystem>("Simple System"); // 往方程系统中添加一个方程，并命名
equation_systems.get_system("Simple System").add_variable("u", FIRST); // 给方程添加变量
equation_systems.add_system<ExplicitSystem>("Complex System"); // 再次添加一个方程
equation_systems.get_system("Complex System").add_variable("c", FIRST); // 添加c变量
equation_systems.get_system("Complex System").add_variable("T", FIRST); // 添加T变量
equation_systems.get_system("Complex System").add_variable("dv", SECOND, MONONIAL); // 添加dv变量
equation_systems.init(); // 初始化方程系统的数据结构
mesh.print_info(); // 在屏幕上打印网格信息
equation_systems.print_info(); // 在屏幕上打印方程系统信息
equation_systems.write(argv[1],WRITE); // 将方程信息输出到文件中
equation_systems.clear(); // 清除方程的数据结构
equaiton_systems.read(argv[1],READ); // 读入刚才的文件, 这三步可用于周期性地存储结果，并接着计算
Intro 3
这个例子演示了怎样求解Possion方程。同样展示了矩阵组装、解析解对比、单元迭代等功能。
头文件：1
2
3
4
5
6
7
8
9
10
11
// 定义有限元刚度矩阵和解向量
源码：1
2
3
4
5
6
7
8
9
10MeshTools::Generation::build_square(mesh, 15,15,-1.0,1.0,-1.0,1.0,QUAD9); // 设定生成网格的计算域、网格密度和类型
EquationSystems equation_systems(mesh);
equation_systems.add_system<LinearImplicitSystem>("Possion");
equation_systems.get_system("Possion").add_variable("u", SECOND);
// 指定一个指向矩阵组装功能的指针，具体的组装方法见函数实现
equation_systems.get_system("Possion").attatch_assemble_function (assemble_possion);
// 求解Possion方程，这将组装线性系统，并调用数值求解器，对于PETSc来说，可以在命令行中指定求解器 -ksp_type cg
equation_systems.get_system("Possion").solve();
具体的组装方法：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52const MeshBase & mesh = es.get_mesh(); // 得到mesh对象的const引用
const unsigned int dim = mesh.mesh_dimension(); // 得到dimension
LinearImplicitsystem & system = es.get_system<LinearImplicitsystem>("Possion"); // 得到方程系统的引用
const DofMap & dof_map = system.get_dof_map(); // 得到自由度映射
FEType fe_type = dof_map.variable_type(0); // 得到系统的第一个且这里是唯一一个变量的有限单元类型的引用
UniquePtr<FEBase> fe (FEBase::build(dim,fe_type)); // 基于有限单元类型创建一个有限元对象
QGauss qrule (dim, FIFTH); // 使用5次高斯积分准则
fe -> attatch_quadrature_rule (&qrule); //告诉有限元对象使用该高斯积分准则
UniquePtr<FEBase> fe_face (FEBase::build(dim,fe_type)); // 创建一个特殊的有限元对象，用于边界积分
QGauss qface (dim-1, FIFTH); // 边界上的积分准则的dimension比单元的小1
fe_face -> attatch_quadrature_rule (&qface); //告诉有限元对象使用该高斯积分准则
const std::vector<Real> & JxW = fe->get_JxW(); // 得到Jacobian矩阵与积分权重的乘积
const std::vector<Point> & q_point = fe -> get_xyz(); // 得到单元积分点的坐标
const std::vector<std::vector<Real> > & phi = fe -> get_phi(); // 得到积分点上的形函数
const std::vector<std::vector<RealGradient> > & dphi = fe -> get_dphi(); // 得到积分点上的形函数梯度
DenseMatrix<Number> Ke; // 定义存储单元刚度矩阵的数据结构
DenseVector<Number> Fe; // 定义存储单元右端项向量的数据结构
std::vector<dof_id_type> dof_indices; // 存储自由度的全局标识
// 下面是对所有单元进行循环，计算单元刚度矩阵和右端项对全局的贡献
MeshBase::const_element_iterator el = mesh.active_local_elements_begin(); // 起始单元
const MeshBase::const_element_iterator end_el = mesh.active_local_element_end(); // 结束单元
for (; el != end_el; ++el) // for循环
{ const Elem * elem = *el; // 存储一下当前单元的指针，便于后续的语法表示
dof_map.dof_indices(elem, dof_indices); // 得到当前单元的全局自由度标识，从而知道该单元往哪个全局自由度上做贡献
fe->reinit(elem); // 计算当前单元的特定信息，包括积分点的位置和形函数及其梯度。
Ke.resize(dof_indices.size(),dof_indices.size()); // 在累加之前先置0
Fe.resize(dof_indices.size()); //在累加之前先置0
// 然后在该单元的所有积分点上进行循环
for (unsigned int qp=0; qp < qrule.n_points(); qp++){
// 下面是构建单元刚度矩阵，包括一个嵌套循环
for (std::size_t i=0; i< phi.size(); i++)
for(std::size_t j=0; j< phi.size(); j++)
{
Ke(i,j)+=JxW[qp]*(dphi[i][qp]*dphi[j][qp]);
}
// 下面是构建单元右端项
for (std::size_t i=0; i < phi.size(); i++)
Fe(i) += JxW[qp]*fxy*phi[i][qp]; //fxy是源项
}
// 还要施加边界条件，该部分详见原代码
system.matrix->add_matrix(Ke, dof_indices); // 将单元刚度矩阵的贡献叠加到全局上
system.rhs->add_vector(Fe, dof_indices); // 将单元右端项的贡献叠加导全局上
}
Intro 4
该例子是基于例子3，展示了通过很小的改动来实现与维度无关的编程。
同时引出PerfLog类来查看性能，从而确定瓶颈所在以实现优化。
头文件：1
2
3
源代码：1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19GetPot command_line (argc, argv); // 创建一个GetPot对象来解析命令行
int dim = 2;
if(command_line.search(1, "-d"))
dim = command_line.next(dim); // 读取命令行中设定的dim
// 根据不同的dimension来生成网格
MeshTools::Generation::build_cube (mesh,
ps,
(dim>1) ? ps : 0,
(dim>2) ? ps : 0,
-1., 1.,
-halfwidth, halfwidth,
-halfheight, halfheight,
(dim==1) ? EDGE2 :
((dim == 2) ? QUAD4 : HEX8));
perf_log.push("elem init"); // 开始记录形函数初始化过程
perf_log.pop("elem init"); // 结束记录形函数初始化过程，push和pop成对出现，否则会报错
Intro 5
第五个例子展示了怎样方便调整积分准则。
源码：1
2QuadratureType quad_type = INVALID_Q_RULE; // 定义积分类型
quad_type = static_cast<QuadratureType>(std::atoi(argv[2])); // 从命令行中读取积分类型
