Gmsh example 1

什么是Gmsh?

Gmsh是一个开源的具有内置CAD引擎和后期处理子的三维有限元网格生成器. 它的目的是提供一个快速的, 轻量级的, 用户友好的网格工具, 并且具有参数化输入和后期可视化能力. 目前Gmsh具有四个模块: 几何模块, 网格模块, 解子模块和后期处理模块, 且所有模块的使用都可以通过两种方式: 命令行和GUI.
有限元网格生成工具有很多种, 比如MeshGen, Hypermesh, 以及一个MATLAB小工具DistMesh, 等等. 可选择性很多, 但很多时候把时间浪费在学习工具的使用上面是不明智的, 工具有一个就够了. 还有一个更实际的问题, 作为一个没有固定收入的学生来说, 一个开源免费实用的软件是多么重要!
Gmsh主页http://gmsh.info/,有Linux, Mac and Windows版本可用, 并且在下载的安装包的tutorial目录包含一个介绍关键概念和一些实例的教程, 此处应该点赞.

第一个例子.

在学习有限元的时候, 我们需要一些标准的网格, 比如在二维区域[0,1]*[0,1]上的结构网格, 来实现有限元的程序. 一般情况下, 结点坐标和单元是一个网格最基本的内容, 也是编写程序必不可少的数据. 很多情况下, 为了验证算法在不同网格下的收敛情况, 生成具有相应的性质的网格又是必须的. 比如下面的SCILAB代码可以生成[0,1]*[0,1]上如图所示的标准结构网格.

gmsh_example_1_structure_mesh

function [ node,Nn, elem,Ne ] = mesh_gen2d( N )
//MESH_GEN - Divide the domain [0,1]*[0,1] into 2*N*N triangles
//   N - Number in every dimension
  h=1.0/N;
  Ne = 2*N*N;       // numerical of node
  Nn=(N+1)*(N+1);   // numerical of elem
  node = zeros(Nn,2); elem = zeros(Ne, 3);
  for i=0:N         // node 
      for j=0:N
          node(i*(N+1)+j+1, 1:2)=[j*h, i*h];
      end
  end
  ie=1;             // elem
  for i=0:N-1
      in=i*(N+1)+1;
      for j=0:N-1
          elem(ie,:) = [in+j, in+j+1, in+j+N+1];
          elem(ie+1,:)=[in+j+1, in+j+N+2, in+j+N+1];
          ie=ie+2;
      end
  end
endfunction

我们现在比较关心的是如何利用Gmsh来生成一个简单的结构网格, 严格来说Gmash生成的网格都非结构的, 因为它无法保证所有网格具有相同的大小和形状. Gmsh是三维网格生成器, 所以三维空间里的点, 线, 面, 体都补看作网格里的单元, 例如可以通过如下语句定义一个点元:

Point(1) = {0.0, 0.0, 0.0, lc2};

如果仅仅是这样的话, 在之后生成的.msh文件中是不会把这个点看作网格中的一个物理单元记录下它的数据的, 即不会在数据块$Elements $EndElements中出现该点信息. 不过我们可以通过如下语句把这个点元升级为一个物理点元:

 Physical Point(1) = {1};

那么数据块$Elements $EndElements中就会多出如下该点(作为一个点元出现)的信息:

$Elements
1 15 2 1 1 1
****
$EndElements

关于单元和物理实体的区别, 以及引入物理实体这个概念的好处在官方手册6.2节http://gmsh.info/doc/texinfo/gmsh.html有比较详细的描述.

这里给出一个例子用来生成一个近似的结构网格, 请看如下gmsh文件ex1.geo:

///////////////////////////////////////////////////////////////////////////////
lc = 0.5;
lc2 = 0.5;  //lc,lc2 - Control size of elements around points
Point(1) = {0.0, 0.0, 0.0, lc2};
Point(2) = { 1.0, 0.0, 0.0, lc};
Point(3) = { 1.0,  1.0, 0.0, lc};
Point(4) = {0.0,  1.0, 0.0, lc};

Line(1) = {1,2};
Line(2) = {2,3};
Line(3) = {3,4};
Line(4) = {4,1};

Line Loop(5) = {1,2,3,4};

Plane Surface(6) = {5};
Physical Surface(1) = {6} ;
///////////////////////////////////////////////////////////////////////////////

因为我们想要的是网格单元的数据, 所以定义了一个物理面surface, 即语句Physical Surface(1) = {6}, 其目的就是只得到区域上结点坐标和三角单元信息. 读者可以自己在Gmsh里Mesh模块中生成一下.msh文件, 看是不是只有三角单元的数据信息. 如果我们想让Point(1)附近的网格比较细, 这在处理一些具有低正则解或奇异解的问题时是合理的, 那么只需要把lc2的值改为0.01或者更小, 就会生成另外一个非结构网格, 为了比较, 笔者附上了两种情况下生成的不同网格图像.

gmsh_example_1_structure_mesh2 gmsh_example_1_structure_mesh3