数字旗手

引子 这个例子的亮点：介绍了最简单的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

本文是对Mathematica有限元分析与工程应用一书的学习笔记。 三角形单元适用于具有复杂和不规则边界形状的问题，是一种常见的网格离散方式。 双线性三角形单元 局部坐标系 之前的杆单元和桁架元的局部坐标系容易建立，即建在其自身上即可。而三角形单元的局部坐标系显然不能这样建立，其常采用一套无量纲的自然坐标系——面积坐标，如下图所示： 三角形123内部有一任意点P，P与顶点1、2、3组成3个子三角形，每个子三角形的面积与总面积之比记为$L_i$，即P点的面积坐标为$(L_1,L_2,L_3)$。 因为三个坐标相加为1，所以仅有两个独立的自然坐标，所以可以简记为： 注意到：面积坐标还有如下

本文是对Mathematica有限元分析与工程应用一书的学习笔记。 桁架元的特点 平面桁架元是既有局部坐标又有整体坐标的二维有限元，因此比起之前的杆单元，需要多一步坐标变换。 桁架元示意图如下： 指定整体坐标系为X-Y，局部坐标系为x-y。则两者之间的转换关系为： 即： 其中： 局部坐标系下的有限元方程为： 为了把有限元方程从局部坐标系变换到整体坐标系，可通过转换矩阵： 所以： 又因为转换矩阵T满足如下关系(可实际计算验证一下)： 所以： 所以整体坐标系的刚度矩阵与局部坐标系的刚度矩阵关系为： 由于局部坐标系下的单元刚度矩阵为： 其中$k=\frac{EA}{L}$。那