引言
用于枚举基函数(形函数)的类称为DoFHandler。
要在网格上定义自由度,先要创建一个有限元对象,然后通过DoFHandler::distribute_dofs 函数传给DoFHandler。(distribute DoFs 代表了枚举基函数的过程。)DoFHandler 会提供总体装配信息。
稀疏矩阵
DoFHandler 默认以随机方式枚举自由度,因此相应的稀疏矩阵也是随机排列的。一般来说随机性不影响求解,但某些预处理例如SSOR在特定排列下会更快些。
程序
网格库
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/manifold_lib.h>
把自由度关联到单元上的库
#include <deal.II/dofs/dof_handler.h>
描述双线性单元(包括1至3维)的库
#include <deal.II/fe/fe_q.h>
处理自由度的若干工具
#include <deal.II/dofs/dof_tools.h>
稀疏矩阵
#include <deal.II/lac/sparse_matrix.h>
动态(即瞬时)稀疏矩阵
#include <deal.II/lac/dynamic_sparsity_pattern.h>
重排列
#include <deal.II/dofs/dof_renumbering.h>
输出
#include <fstream>
命名空间
using namespace dealii;
网格生成
void make_grid(Triangulation<2> &triangulation)
{
const Point<2> center(1, 0);
const double inner_radius = 0.5, outer_radius = 1.0;
GridGenerator::hyper_shell(
triangulation, center, inner_radius, outer_radius, 5);
for (unsigned int step = 0; step < 3; ++step)
{
for (auto &cell : triangulation.active_cell_iterators())
for (unsigned int v = 0; v < GeometryInfo<2>::vertices_per_cell; ++v)
{
const double distance_from_center =
center.distance(cell->vertex(v));
if (std::fabs(distance_from_center - inner_radius) < 1e-10)
{
cell->set_refine_flag();
break;
}
}
triangulation.execute_coarsening_and_refinement();
}
}
创建DoFHandler。首先要创建一个拉格朗日单元FE_Q,它的唯一参数表示多项式次数,双线性为1。接着把它传给DoFHandler。
void distribute_dofs(DoFHandler<2> &dof_handler)
{
const FE_Q<2> finite_element(1);
dof_handler.distribute_dofs(finite_element);
然后考虑稀疏矩阵。存储非零元素位置的类是SparsityPattern。这个类有一个毛病,必须先给定单行最大存储数。二维时可以用DoFHandler::max_couplings_between_dofs()
函数自动获取,但三维时总是会过分估计。为此我们可以用一个DynamicSparsityPattern,再复制到SparsityPattern里。DynamicSparsityPattern 需要输入矩阵大小。
DynamicSparsityPattern dynamic_sparsity_pattern(dof_handler.n_dofs(),
dof_handler.n_dofs());
用自由度填充DynamicSparsityPattern。
DoFTools::make_sparsity_pattern(dof_handler, dynamic_sparsity_pattern);
最后曲线救国,把位置信息复制到SparsityPattern。
SparsityPattern sparsity_pattern;
sparsity_pattern.copy_from(dynamic_sparsity_pattern);
把SparsityPattern输出到图片。
std::ofstream out("sparsity_pattern1.svg");
sparsity_pattern.print_svg(out);