数字旗手

电气化、自动化、数字化、智能化、智慧化

0%

引子 本例主要讲解MeshWorker框架和间断Galerkin方法,即DG。 * 用间断Galerkin法离散线性对流方程 * 使用MeshWorker::loop()来组装系统矩阵 本例主要关心的就是间断Galerkin法的循环,这相当复杂,因为必须分辨边界、常规内部边和有悬点的内部边。MeshWorker框架能够对所有的单元和边进行标准循环,它将分辨过程隐藏在了内部。 使用MeshWorker需要手动做两件事:一是针对特定问题写内部积分器,二是从该命名空间中选择类,然后将它们组合起来来解决问题。 要求解的问题是线性对流方程: 边界条件是: 这是入流边界,定义是: 这个方程
Read more »

引子 这个小例子主要演示更高阶的映射。映射的意思是参考单元和真实单元之间的变换,之前的例子中默认使用线性或双线性映射。这里的线性映射不是形函数用线性近似,同时也不是Cn单元的概念,Cn单元代表解的最高阶导数的概念。如果边界是弯曲的话,用线性近似就可能不充分。如果使用分段二次抛物线来近似的话,就说这是二次或Q2近似。如果使用分段三次多项式近似的话,就称为Q3近似。依次类推。 这里是用逼近圆周率来说明映射问题。用两种方法:一种是对半径为1的圆进行三角剖分,然后积分,如果完全是圆,那么它的面积精确为圆周率,但因为是多项式近似,所以肯定不精确,这里展示不同映射下随着网格细化的收敛性。第二种不是计算面积
Read more »

引子 这个例子的亮点:介绍了最简单的Laplace方程的组装矩阵和右端项的一键生成函数。 这个例子来求解仅考虑Neumann边界条件的Laplace问题: 如果有解的话,解必须满足协调条件: 这里考虑计算域是以原点为圆心,半径为1的圆,$f=-2,g=1$是满足协调条件的。 虽然协调条件允许有解,但解仍然是不定的:在方程中仅仅解的导数固定,解可以是任意常数。所以需要施加另外一个条件来固定这个常数。可以有以下可能的方法: * 固定离散后的某个节点值为0或其他固定值。这意味着一个额外的条件$u_h(x_0)=0$。尽管这是通常的做法,但不是一个好方法,因为我们知道Laplace方程的解是在
Read more »

引子 本例将要完成以下目标: * 求解对流方程$\beta \cdot \nabla u = f$ * 使用多线程求解 * 设计一个简单的细化准则 方程离散 对流方程中$\beta$是一个描述对流方向和速度的矢量场,可能与空间位置相关,$f$是源项,$u$是解。这个方程求解的物理过程可以是浓度u在给定流场下的传输,没有扩散过程,但有源。 明显地,在入流侧,上述方程需要在边界上获得补充: 其中,$\Omega_-$代表边界的入流部分,其定义为: 其中,${\mathbf n}({\mathbf x})$是点$x$上的外法线。这种定义方式是非常直观的:因为如果外法线是朝外的,那么在入
Read more »

引子 真实生活中,大多数的偏微分方程都是一组方程,相应地,解也通常是矢量场。跟之前单方程的求解以及解是标量场相比,这种问题只是在组装矩阵和右端项时有些复杂,其他都一样。 本例中求解的是弹性问题: 其中,$c_{ijkl}$是刚度系数,通常与空间坐标相关。 弹性方程是Laplace方程的一种扩散形式,其解$u_l$是矢量场,表示一个刚体在力的作用下的位移。力$f_i$也是矢量场。 $c_{ijkl}$是一个四阶张量,共有$3^4=81$个分量,但其实只有21个分量是独立的,这是因为: 首先因为$\sigma_{ij}$和$\epsilon_{ij}$都是对称张量,导致$c_{ijkl}=c_{
Read more »

引子 在本例中,将会着眼于以下两方面: 1. 验证程序的正确性,生成收敛性统计表格; 2. 对于Helmholtz方程施加非齐次Neumann边界条件。 另外还有一些小的优化点。 验证程序正确性 也许从来不会有任何一个有限元程序一开始就是正确的,所以找到方法来验证计算的解是否正确就很有必要。通常选择已知精确解析解,并且比较精确解析解和计算离散解两者之间差别来求证。如果随着误差次数提高,两者之间差别逐渐趋于0,就说明程序的正确性。deal.II中就提供了这样一个函数:VectorTools::integrate_difference(),它提供了很多种范数的计算: 这些公式也适用于矢
Read more »

引子 本例主要着眼于处理局部细化的网格。如果临近单元细化多次以后,单元界面上的格点可能在另一边就不平衡,称为“悬点”。为了保证全局解在这些格点上也是连续的,必须对这些节点上的解的值施加一些限制,相应的核心的类是ConstraintMatrix。 程序解析 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 #include #include #include
Read more »

引子 此例没有介绍革命性的功能,但有很多对前面例子的“微创新”,包括: * 在不断细化的网格上的计算。数值计算通常要在不同的网格上进行,这样才能感受到精度。而且deal.II支持自适应网格,虽然这个例子中没有用到,但基础在这 * 读入非规则网格数据 * 计算优化 * debug模式,使用Assert宏 * 变系数Possion方程,使用预条件迭代求解器 这里要求解的方程是: 如果$a(x)$是常系数,那么就成了Possion方程,如果它是空间相关的变系数,方程就复杂一些了。 还是得先写出方程的弱形式: 程序解析 以下是头文件们: 1 2 3 4 5 6 7 8 9 10 11
Read more »

引子 deal.II有一个特性,叫作”维度无关的编程“。前面的例子中,很多类都有一个尖括号括起的数字的后缀。 这意味着该类能作用在不同的维度上,而不同维度的计算代码基本相同,这能显著地减少重复代码。这正是C++的模板template的拿手好戏。 在Step4中,将显示程序怎样维度无关的编程,具体是将step3中的Laplace问题同时在二维和三维上求解,以及右端项不再是常量、边界值不再为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 #include
Read more »

引子 这是使用有限元法进行具体计算的第一个算例,求解的是一个简化的Possion方程,其在边界上为0,而右端项不为0,即: 求解域是单位正方形$\Omega=[0,1]^2$,其上的网格划分在step1和step2中已经涉及。 这里也仅仅计算特例$f(x)=1$,更一般的情形详见step4。 推导 首先需要得到上述方程的弱形式,即在方程两侧左乘一个测试函数$\phi$并在计算域上积分(注意是左乘,而不是右乘): 然后利用高斯散度公式进行分部积分,可得: 同时,因为测试函数必须满足相同的边界条件,即在边界上$\phi=0$,所以上式变为: 这里应用了通常的简写形式:$(a,b)=\in
Read more »