数字旗手

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

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

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

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

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

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

引子 本例主要着眼于处理局部细化的网格。如果临近单元细化多次以后，单元界面上的格点可能在另一边就不平衡，称为“悬点”。为了保证全局解在这些格点上也是连续的，必须对这些节点上的解的值施加一些限制，相应的核心的类是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

引子 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