求解偏微分方程开源有限元软件deal.II学习--Step 10
求解偏微分方程开源有限元软件deal.II学习--Step 10
引子
这个小例子主要演示更高阶的映射。映射的意思是参考单元和真实单元之间的变换,之前的例子中默认使用线性或双线性映射。这里的线性映射不是形函数用线性近似,同时也不是Cn单元的概念,Cn单元代表解的最高阶导数的概念。如果边界是弯曲的话,用线性近似就可能不充分。如果使用分段二次抛物线来近似的话,就说这是二次或Q2近似。如果使用分段三次多项式近似的话,就称为Q3近似。依次类推。
这里是用逼近圆周率来说明映射问题。用两种方法:一种是对半径为1的圆进行三角剖分,然后积分,如果完全是圆,那么它的面积精确为圆周率,但因为是多项式近似,所以肯定不精确,这里展示不同映射下随着网格细化的收敛性。第二种不是计算面积,而是计算周长,那么圆周率就大约是周长的一半。
程序解析
1 2 3 4 5 6 7 8 9 10 11 12 |
#include <deal.II/base/quadrature_lib.h> #include <deal.II/base/convergence_table.h> #include <deal.II/grid/grid_generator.h> #include <deal.II/grid/tria_accessor.h> #include <deal.II/grid/tria_iterator.h> #include <deal.II/grid/manifold_lib.h> #include <deal.II/grid/tria.h> #include <deal.II/grid/grid_out.h> #include <deal.II/dofs/dof_handler.h> #include <deal.II/dofs/dof_accessor.h> #include <deal.II/fe/fe_q.h> #include <deal.II/fe/fe_values.h> |
这是之前用过的头文件。
1 |
#include <deal.II/fe/mapping_q.h> |
这是新的头文件。用来声明MappingQ类,进行任意阶数的多项式映射。
1 2 3 |
#include <iostream> #include <fstream> #include <cmath> |
C++的标准头文件。
然后开始定义这个问题的命名空间。
1 2 3 |
namespace Step10 { using namespace dealii; |
先定义一种非常精确的圆周率值用于后面比较:
1 |
const long double pi = 3.141592653589793238462643; |
下面是输出。这里因为是个小程序,没有用类,而是用函数模板,参数是空间维度。
1 2 3 4 5 6 7 8 9 10 |
template <int dim> void gnuplot_output() { std::cout << "Output of grids into gnuplot files:" << std::endl << "===================================" << std::endl; Triangulation<dim> triangulation; GridGenerator::hyper_ball (triangulation); static const SphericalManifold<dim> boundary; triangulation.set_all_manifold_ids_on_boundary(0); triangulation.set_manifold (0, boundary); |
首先产生一个粗的三角剖分,施加一个合适的边界描述。
1 2 3 4 5 6 7 8 9 10 |
for (unsigned int refinement=0; refinement<2; ++refinement, triangulation.refine_global(1)) { std::cout << "Refinement level: " << refinement << std::endl; std::string filename_base = "ball"; filename_base += '0'+refinement; for (unsigned int degree=1; degree<4; ++degree) { std::cout << "Degree = " << degree << std::endl; const MappingQ<dim> mapping (degree); |
然后使用不同的映射,分别是Q1、Q2和Q3。当分段线性映射,即Q1时,可以直接使用MappingQ1这个类。Attention!!MappingQ1也是很多函数和类默认使用的映射方式,如果不明确指定映射方式的话。
1 2 3 4 5 6 7 8 9 10 11 12 |
GridOut grid_out; GridOutFlags::Gnuplot gnuplot_flags(false, 30); grid_out.set_flags(gnuplot_flags); std::string filename = filename_base+"_mapping_q"; filename += ('0'+degree); filename += ".dat"; std::ofstream gnuplot_file (filename.c_str()); grid_out.write_gnuplot (triangulation, gnuplot_file, &mapping); } std::cout << std::endl; } } |
然后就是使用Gnuplot进行输出。
下面就进入该算例的主要部分:圆周率的计算。
因为这里的圆的半径是1,所以圆的面积的计算就是对常数1在整个计算域上积分:
计算开始:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 |
template <int dim> void compute_pi_by_area () { std::cout << "Computation of Pi by the area:" << std::endl << "==============================" << std::endl; const QGauss<dim> quadrature(4); for (unsigned int degree=1; degree<5; ++degree) { std::cout << "Degree = " << degree << std::endl; Triangulation<dim> triangulation; GridGenerator::hyper_ball (triangulation); static const SphericalManifold<dim> boundary; triangulation.set_all_manifold_ids_on_boundary (0); triangulation.set_manifold(0, boundary); const MappingQ<dim> mapping (degree); |
这里选择的积分精度足够大,保证在这里任意映射下都能正确求解。
1 2 |
const FE_Q<dim> dummy_fe (1); DoFHandler<dim> dof_handler (triangulation); |
这里创建了虚假的有限单元和自由度句柄,因为这里不需要知道在积分点上有限元上的值,只是为了知道那个积分权重。注意,这里有限单元的形函数次数是1,再次说明线性形函数跟线性映射不是一个概念。
1 2 |
FEValues<dim> fe_values (mapping, dummy_fe, quadrature, update_JxW_values); |
这里传递给FEValues的第一个参数是mapping对象,之前的例子中这个参数是省略的,默认使用MappingQ1这样的对象。
1 |
ConvergenceTable table; |
这里再创建一个收敛性表格,来看收敛速率。
然后开始对细化步数循环:
1 2 3 4 5 6 |
for (unsigned int refinement=0; refinement<6; ++refinement, triangulation.refine_global (1)) { table.add_value("cells", triangulation.n_active_cells()); dof_handler.distribute_dofs (dummy_fe); long double area = 0; |
然后对所有单元进行循环:
1 2 3 4 5 6 7 8 9 10 11 12 |
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); for (; cell!=endc; ++cell) { fe_values.reinit (cell); for (unsigned int i=0; i<fe_values.n_quadrature_points; ++i) area += fe_values.JxW (i); } table.add_value("eval.pi", static_cast<double> (area)); table.add_value("error", static_cast<double> (std::fabs(area-pi))); } |
并且存储进table中。
然后开始计算收敛性:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 |
table.omit_column_from_convergence_rate_evaluation("cells"); table.omit_column_from_convergence_rate_evaluation("eval.pi"); table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2); table.set_precision("eval.pi", 16); table.set_scientific("error", true); table.write_text(std::cout); std::cout << std::endl; } } 下面的计算方法不是计算面积,而是计算周长: ```cpp template <int dim> void compute_pi_by_perimeter () { std::cout << "Computation of Pi by the perimeter:" << std::endl << "===================================" << std::endl; const QGauss<dim-1> quadrature(4); |
注意,上面的QGauss的维度是dim-1,这是因为是对边积分,而不是对体单元。
1 2 3 4 5 6 7 8 9 10 11 |
for (unsigned int degree=1; degree<5; ++degree) { std::cout << "Degree = " << degree << std::endl; Triangulation<dim> triangulation; GridGenerator::hyper_ball (triangulation); static const SphericalManifold<dim> boundary; triangulation.set_all_manifold_ids_on_boundary (0); triangulation.set_manifold (0, boundary); const MappingQ<dim> mapping (degree); const FE_Q<dim> fe (1); DoFHandler<dim> dof_handler (triangulation); |
这些跟上面的一样。
1 2 3 4 5 6 7 8 |
FEFaceValues<dim> fe_face_values (mapping, fe, quadrature, update_JxW_values); ConvergenceTable table; for (unsigned int refinement=0; refinement<6; ++refinement, triangulation.refine_global (1)) { table.add_value("cells", triangulation.n_active_cells()); dof_handler.distribute_dofs (fe); |
Attention!这里是创建的FEFaceValues而不是FEValues对象。
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 |
typename DoFHandler<dim>::active_cell_iterator cell = dof_handler.begin_active(), endc = dof_handler.end(); long double perimeter = 0; for (; cell!=endc; ++cell) for (unsigned int face_no=0; face_no<GeometryInfo<dim>::faces_per_cell; ++face_no) if (cell->face(face_no)->at_boundary()) { fe_face_values.reinit (cell, face_no); for (unsigned int i=0; i<fe_face_values.n_quadrature_points; ++i) perimeter += fe_face_values.JxW (i); } table.add_value("eval.pi", static_cast<double> (perimeter/2.)); table.add_value("error", static_cast<double> (std::fabs(perimeter/2.-pi))); } table.omit_column_from_convergence_rate_evaluation("cells"); table.omit_column_from_convergence_rate_evaluation("eval.pi"); table.evaluate_all_convergence_rates(ConvergenceTable::reduction_rate_log2); table.set_precision("eval.pi", 16); table.set_scientific("error", true); table.write_text(std::cout); std::cout << std::endl; } } } |
Attention!!这里积分时,只有边界上的才能被叠加进去,所以需要有一个判断是否边界的语句。
然后是main函数:
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 |
int main () { try { std::cout.precision (16); Step10::gnuplot_output<2>(); Step10::compute_pi_by_area<2> (); Step10::compute_pi_by_perimeter<2> (); } catch (std::exception &exc) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Exception on processing: " << std::endl << exc.what() << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } catch (...) { std::cerr << std::endl << std::endl << "----------------------------------------------------" << std::endl; std::cerr << "Unknown exception!" << std::endl << "Aborting!" << std::endl << "----------------------------------------------------" << std::endl; return 1; } return 0; } |
计算结果
上面展示的是粗网格上的结果。实线是网格,虚线是精确的圆。
可以看出,即使在粗网格上,二次和三次映射都取得了很好的吻合结果。
然后是计算精度:
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 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 |
Computation of Pi by the area: ============================== Degree = 1 cells eval.pi error 5 1.9999999999999993 1.1416e+00 - 20 2.8284271247461894 3.1317e-01 1.87 80 3.0614674589207178 8.0125e-02 1.97 320 3.1214451522580520 2.0148e-02 1.99 1280 3.1365484905459393 5.0442e-03 2.00 5120 3.1403311569547534 1.2615e-03 2.00 Degree = 2 cells eval.pi error 5 3.1045694996615865 3.7023e-02 - 20 3.1391475703122271 2.4451e-03 3.92 80 3.1414377167038303 1.5494e-04 3.98 320 3.1415829366419015 9.7169e-06 4.00 1280 3.1415920457576911 6.0783e-07 4.00 5120 3.1415926155921139 3.7998e-08 4.00 Degree = 3 cells eval.pi error 5 3.1410033851499310 5.8927e-04 - 20 3.1415830393583861 9.6142e-06 5.94 80 3.1415925017363837 1.5185e-07 5.98 320 3.1415926512106722 2.3791e-09 6.00 1280 3.1415926535525962 3.7197e-11 6.00 5120 3.1415926535892140 5.7923e-13 6.00 Degree = 4 cells eval.pi error 5 3.1415871927401127 5.4608e-06 - 20 3.1415926314742437 2.2116e-08 7.95 80 3.1415926535026228 8.7170e-11 7.99 320 3.1415926535894529 3.4036e-13 8.00 1280 3.1415926535897927 2.9187e-16 10.19 5120 3.1415926535897944 1.3509e-15 -2.21 Computation of Pi by the perimeter: =================================== Degree = 1 cells eval.pi error 5 2.8284271247461898 3.1317e-01 - 20 3.0614674589207178 8.0125e-02 1.97 80 3.1214451522580520 2.0148e-02 1.99 320 3.1365484905459393 5.0442e-03 2.00 1280 3.1403311569547525 1.2615e-03 2.00 5120 3.1412772509327729 3.1540e-04 2.00 Degree = 2 cells eval.pi error 5 3.1248930668550594 1.6700e-02 - 20 3.1404050605605449 1.1876e-03 3.81 80 3.1415157631807014 7.6890e-05 3.95 320 3.1415878042798617 4.8493e-06 3.99 1280 3.1415923498174534 3.0377e-07 4.00 5120 3.1415926345932004 1.8997e-08 4.00 Degree = 3 cells eval.pi error 5 3.1414940401456057 9.8613e-05 - 20 3.1415913432549156 1.3103e-06 6.23 80 3.1415926341726914 1.9417e-08 6.08 320 3.1415926532906893 2.9910e-10 6.02 1280 3.1415926535851360 4.6571e-12 6.01 5120 3.1415926535897203 7.2845e-14 6.00 Degree = 4 cells eval.pi error 5 3.1415921029432576 5.5065e-07 - 20 3.1415926513737600 2.2160e-09 7.96 80 3.1415926535810712 8.7218e-12 7.99 320 3.1415926535897594 3.3668e-14 8.02 1280 3.1415926535897922 1.0617e-15 4.99 5120 3.1415926535897931 1.0061e-16 3.40 |
可以看出,随着映射次数提高,计算精度也在提高,而Q1映射在最细网格上的精度还不如Q3映射在最粗网格上的精度。最后一栏是收敛阶,可以看出阶数达到了h2p
。
求解偏微分方程开源有限元软件deal.II学习--Step 10相关推荐
- 求解二阶偏微分方程c语言,科学网—求解偏微分方程开源有限元软件deal.II学习--Step 3 - 亓欣波的博文...
完整版见:qixinbo.info. deal.II的程序结构 deal.II采用面向对象编程,其中包含了很多的Modules,各自实现不同的功能,并有机地结合起来.如上图所示.具体为:Triangu ...
- 流固耦合开源软件precice安装笔记(包括开源CFD软件OpenFOAM、插件swak4Foam,开源有限元软件CalculiX、deal.II)
安装环境:Ubuntu 20.04 LTS 1. 安装Python发行版Anaconda 可在Anaconda官网下载安装包,下载完成后在下载目录中鼠标右键打开终端,键入: bash Anaconda ...
- 开源的python有限元软件_想从Abaqus转用开源有限元软件有什么好的推荐吗?
推荐FreeCAD,Code-Aster,Calculix. (1)FreeCAD据说华为在尝试于FreeCAD基础上进行CAE开发,FreeCAD本是一款开源的CAD软件,具备一些简单的有限元分析功 ...
- 编译开源有限元软件calculix
calculix的编译过程 第一步:下载安装cygwin,选择gcc-core gcc-g++ gfortran perl注意版本一致,且为最低版本 第二步:修改makefile文件,因为最后没有生成 ...
- python开源流程图软件_适用于Linux的10种最佳流程图和图表软件
图表是我们联系信息并处理其重要性的好方法. 它们有助于沟通关系和抽象信息,并使我们可视化概念. 流程图和图表工具可用于从基本工作流程图到复杂网络图 ,组织图, BPMN ( 业务流程模型和表示法 ), ...
- 【哈工大软件构造】学习笔记10 第十章、第十一章、第十二章
目录 第十章 面向可维护性的构造技术 1 软件维护和演化 2 可维护性的度量 3 模块化设计和模块性准则 模块划分的五个准则 模块设计的五个原则 耦合度和聚合度 4 OO设计准则:SOLID SRP ...
- 开源的python有限元软件_python有限元
广告关闭 腾讯云11.11云上盛惠 ,精选热门产品助力上云,云服务器首年88元起,买的越多返的越多,最高返5000元! 节点坐标和单元的拓扑信息暂由有限元前处理软件导出后由python读入#nodes ...
- 【机器人学】机器人开源项目KDL源码学习:(5)KDL如何求解几何雅克比矩阵
这篇文章试图说清楚两件事:1. 几何雅克比矩阵的本质:2. KDL如何求解机械臂的几何雅克比矩阵. 一.几何雅克比矩阵的本质 机械臂的关节空间的速度可以映射到执行器末端在操作空间的速度,这种映射可以通 ...
- 免费开源Scada软件 RapidScada学习记录
免费开源Scada软件 RapidScada学习记录 由于工作需要,需要采集污水处理厂的数据,比如溶解氧,ph等数据,并且要求远程控制功能,设备数据异常手机推送功能,能通过网页查看摄像头当前设备及工作 ...
- 深度学习(Deep Ritz,Galerkin,PINN)求解偏微分方程(PDE)解读
深度学习在求解偏微分方程(PDE)中有三种主要方法: Deep Ritz, Galerkin 和 PINN. Deep Ritz 方法是通过使用深度神经网络来逼近解析解来求解 PDE 的. 它使用了 ...
最新文章
- Android开发历程_18(XML文件解析)
- Linux之系统文件管理
- Gmail技巧之无限别名
- Java中的Google协议缓冲区
- AutoMapper用法一瞥
- Spring-web源码解析之Filter-AbstractRequestLoggingFilter
- Python中的一些特殊函数
- 深度学习 《梯度消失和梯度爆炸》
- 利用MMCM IP核产生用户时钟
- linux crontab 详解
- 使用apache benchmark(ab) 测试报错: apr_socket_recv: Connection timed out (110)
- 修改android的avd路径方法
- serialVersionUID作用
- 视频工具mencoder
- idea java maven 打包,idea maven项目 基于idea自己打包方式 以及使用maven插件打包的三种方式...
- 如何用数据进行产品运营
- 系列篇|一文尽览事件相机原理
- CTF—命令执行总结
- 国内云市场,腾讯云、阿里云、华为云,谁能更胜一筹呢?
- 抖音小视频、千图网图片等多平台的微航去水印微信小程序工具解析