数字旗手

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

本文是对Mathematica有限元分析与工程应用一书的学习笔记。 这里的杆单元是总体坐标与局部坐标一致的一维有限元，因此不需要坐标变换，更易分析。 线性杆单元的形函数 线性杆单元的位移函数是坐标x的一次函数，可用以下代码求解： 1 2 3 4 Clear[u1, u2]; u = a1 + a2 x; solut = Solve[{u1 == (u /. x -> 0), u2 == (u /. x -> 1)}, {a1, a2}]; u = u /. solut[[1]] 其中u1和u2是杆端点x=0和x=1处的位移，结果为： 1 u1 + (-u1 + u2) x

引子 在step1中创建了网格，下面就是在网格上定义自由度。此例中使用一阶线性有限元，其自由度的个数与网格的顶点数相关。后面的例子将展示更高次的单元，其上面的自由度与顶点、边、面及cell都有关。 自由度可以理解为形函数中的系数个数，因为它们是未知的，所以称之为未知量或自由度。 定义网格上的自由度很简单，因为deal.II已经内置该功能了，唯一要做的是创建有限元空间。 头文件 1 #include 该头文件将自由度与顶点、线、cell联系起来。 1 #include 该头文件包

Read more »

简介 deal.II是一款开源的求解偏微分方程的有限元软件，它有如下几个特点： * 使用C++编写 * 有多种单元类型 * 可以大规模并行 * 可以自适应网格 * 文档和范例齐全 * 与其他库有良好的接口 安装 deal.II最新版本为8.4.1,可从官网上下载源码，解压后进入源文件目录安装： 1 2 3 4 5 mkdir build cd build cmake -DCMAKE_INSTALL_PREFIX=/path/to/install/dir ../deal.II make install make test 如果期间需要其他依赖如cmake、doxygen

阅读文献时常常会遇到只给图像却没有原始数据的情形，此时如果想要提取数据，就得借助相应软件，目测用的最多的就是GetData，这个是个商业软件，还有个好用的基于web的开源软件WebPlotDigitizer。这里我们基于Mathematica写一套能用于提取图片中数据点的代码。 参考文献(基本思路参考SE上的这个问题，但具体取点和去点方式不同)： Recovering data points from an image Attention：这里的版本强烈依赖于Mathematica的版本，此处使用的是10.4版本，目测应该使用10以上版本，因为低版本中不会出现工具提示条。 实际坐标系与图像