数字旗手

本文来自亓欣波的博客qixinbo.info。 转载请保留上面信息，谢谢！ 引子 本例亮点： * 精确寻找某个节点并通过要求提供该节点上的自由度指标来取出该点上的值。 * 使用线程并行，否则只能串行等待上一步完成。线程是计算的最小单元。 在本例中，将不太关心描述使用deal.II的新方法，而是着眼于书写模块化及可扩展的有限元程序的方法。 这主要是为了考虑先进研究软件的大小和灵活性：使用了先进误差估计概念和自适应解的应用通常相当庞大。而大家的共识是这样的：庞大的程序，如果没有分成更小更独立的组件，将会很快消亡，因为甚至作者也会很快忘记程序内部各个部分之间的依赖关系。数据封装，比如面向对

引子 本例主要讲解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 »