gmsh手册

pip install --upgrade gmsh

// Gmsh GEO教程1
// 几何基础, 基本实体, 物理组

// Gmsh脚本语言中最简单的构造是" 赋值" . 以下命令定义一个新变量" lc" :
lc = 1e-2;
// 这个变量可用于定义Gmsh最简单的" 基本实体" --" 点" . 一个点由一个标签(一个严格为正的整数; 这里是" 1" )唯一标识, 并由四个数字定义:
// 三个坐标(X, Y和Z)以及该点附近的目标网格尺寸(lc):
Point(1) = {0, 0, 0, lc};
// 网格单元尺寸的分布将通过在整个几何中对这些网格尺寸进行插值获得. 另一种指定网格尺寸的方法是使用通用的网格尺寸场(见t10.geo).
// 一种特殊情况是使用背景网格(见t7.geo).
// 如果没有提供目标网格尺寸, 将基于模型的整体大小为模型使用默认的均匀粗略尺寸.

// 我们可以再定义一些点. 所有点都应具有不同的标签:
Point(2) = {.1, 0, 0, lc};
Point(3) = {.1, .3, 0, lc};
Point(4) = {0, .3, 0, lc};

// 曲线是Gmsh的第二种基本实体, 在曲线中, 直线是最简单的. 一条直线由一个标签标识, 并由两个点的标签列表定义. 例如, 在下面的命令中,
// 直线1从点1开始, 到点2结束.
// 注意, 曲线标签与点标签是分开的--因此我们可以将标签" 1" 重用于第一条曲线. 一般来说, Gmsh中的基本实体标签在每个几何维度中必须是唯一的.
Line(1) = {1, 2};
Line(2) = {3, 2};
Line(3) = {3, 4};
Line(4) = {4, 1};

// 第三种基本实体是面. 为了从上面定义的四条曲线中定义一个简单的矩形面, 首先必须定义一个曲线环. 曲线环也由一个标签(在曲线环中唯一)标识,
// 并由一系列相连的曲线的有序列表定义, 每条曲线都关联一个符号(取决于曲线形成环的方向):
Curve Loop(1) = {4, 1, -2, 3};

// 然后我们可以将面定义为曲线环的列表(这里只有一个, 代表外部轮廓, 因为没有孔洞--见t4.geo中的带孔洞面示例):
Plane Surface(1) = {1};

// 在这一步, Gmsh已经知道如何显示矩形面1并对其进行网格划分. 如果我们想将基本几何实体组合成更有意义的组(例如, 定义一些数学属性(" 域" , " 带有诺依曼条件的边界" ),
// 功能属性(" 左机翼" , " 机身" )或材料属性(" 钢" , " 碳纤维" )), 则需要一个可选步骤.

// 在Gmsh中, 这种组被称为" 物理组" . 默认情况下, 如果定义了物理组, Gmsh在输出文件中只会导出至少属于一个物理组的网格单元.
// (要强制Gmsh保存所有单元, 无论它们是否属于物理组, 设置Mesh.SaveAll=1; 或者在命令行中指定-save_all. )
// 物理组也由标签标识, 即严格为正的整数, 每个维度(0D, 1D, 2D或3D)中的标签必须唯一. 物理组也可以被赋予名称.

// 这里我们定义一个物理曲线, 将左侧, 底部和右侧的曲线组合到一个组中(指定标签5); 以及一个名为" 我的面" 的物理面(带有自动标签), 包含几何面1:
Physical Curve(5) = {1, 2, 4};
Physical Surface("My surface") = {1};

// 现在几何已经完成, 你可以
// - 要么用Gmsh打开这个文件, 在" 网格" 模块中选择" 2D" 来创建网格; 然后选择" 保存" 将其以默认格式保存到磁盘(或者使用" 文件->导出" 导出为其他格式);
// - 要么运行gmsh t1.geo -2在命令行中以批处理模式进行网格划分.

// 你也可以取消注释本脚本中的以下行:
//
// Mesh 2;
// Save "t1.msh";
//
// 这将使Gmsh在每次解析文件时都进行网格划分并保存网格. (要从命令行简单解析文件, 可以使用gmsh t1.geo -)

// 默认情况下, Gmsh将网格保存为最新版本的Gmsh网格文件格式(" MSH" 格式). 你可以通过在图形用户界面, 命令行或脚本中指定具有不同扩展名的文件名,
// 将网格保存为其他网格格式. 例如
//
// // Save "t1.unv";
// 将把网格保存为UNV格式. 你也可以将网格保存为旧版本的MSH格式:
//
// - 在图形用户界面中: 打开" 文件->导出" , 输入你的filename.msh, 然后在下拉菜单中选择版本.
// - 在命令行中: 使用-format选项(例如gmsh file.geo -format msh2 -2).
// - 在.geo脚本中: 添加Mesh.MshFileVersion = x.y; 对于任何版本号x.y.
// - 另一种方法是, 你也可以不明确指定格式, 只需选择带有.msh2或.msh4扩展名的文件名.

// 注意, 从Gmsh 3.0开始, 模型可以使用除默认内置内核之外的其他几何内核构建. 通过指定
//
// SetFactory("OpenCASCADE");
//
// .geo文件中的所有后续命令将由OpenCASCADE几何内核而不是内置内核处理. 不同的几何内核具有不同的功能.
// 使用OpenCASCADE时, 无需通过依次定义4个点, 4条曲线和1个曲线环来定义矩形面, 而是可以直接定义矩形面:
//
// Rectangle(2) = {.2, 0, 0, .1, .3};
//
// 可以使用Boundary或CombinedBoundary运算符访问基础曲线和点.
//
// 参见例如t16.geo, t18.geo, t19.geo或t20.geo获取基于OpenCASCADE的完整示例, 以及examples/boolean获取更多示例.

// Gmsh GEO教程2
// 变换, 拉伸几何, 体积

// 我们首先包含之前的教程文件, 以便将其作为本教程的基础. 包含一个文件相当于将其内容复制粘贴过来:
Include "t1.geo";

// 我们可以用与在t1.geo中相同的方式添加新的点和曲线:
Point(5) = {0, .4, 0, lc};
Line(5) = {4, 5};

// Gmsh还提供了变换(平移, 旋转等)基本实体或基本实体副本的工具. 例如, 点5可以向左移动0.02:
Translate {-0.02, 0, 0} { Point{5}; }

// 还可以围绕(0, 0.3, 0)将其进一步旋转-Pi/4(沿z轴旋转):
Rotate {{0,0,1}, {0,0.3,0}, -Pi/4} { Point{5}; }

// 注意, Gmsh中没有单位: 坐标只是数字--由用户赋予它们含义.

// 点3可以被复制并沿y轴平移0.05:
Translate {0, 0.05, 0} { Duplicata{ Point{3}; } }

// 这个命令创建了一个带有自动分配标签的新点. 可以通过将鼠标悬停在该点上, 在图形用户界面中获取这个标签: 在这种情况下, 新点的标签是" 6" .
Line(7) = {3, 6};
Line(8) = {6, 5};
Curve Loop(10) = {5,-8,-7,3};
Plane Surface(11) = {10};

// 为了自动化工作流程, 不必使用图形用户界面获取新创建实体的标签, 可以直接使用变换命令的返回值. 例如, Translate命令
// 返回一个包含被变换实体标签的列表. 让我们使用以下命令将两个面1和11的副本向右平移:
my_new_surfs[] = Translate {0.12, 0, 0} { Duplicata{ Surface{1, 11}; } };

// my_new_surfs[](注意方括号和命令末尾的; )表示一个列表, 其中包含两个新面的标签
// (查看" 工具->消息控制台" 查看消息):
Printf("新面 '%g' 和 '%g'", my_new_surfs[0], my_new_surfs[1]);

// 在Gmsh中, 列表使用方括号定义(mylist[] = {1, 2, 3};)以及访问其元素(myotherlist[] = {mylist[0], mylist[2]}; mythirdlist[] = myotherlist[];),
// 列表索引从0开始. 要获取列表的大小, 请使用#号: len = #mylist[].

// 注意, 也可以使用圆括号代替方括号, 因此我们也可以写myfourthlist() = {mylist(0), mylist(1)};.

// 体积是Gmsh中的第四种基本实体. 与通过定义曲线环来构建面的方式相同, 必须定义面环(即" 壳" )来构建体积.
// 以下体积没有孔洞, 因此仅由一个面环组成:
Point(100) = {0., 0.3, 0.12, lc};
Point(101) = {0.1, 0.3, 0.12, lc};
Point(102) = {0.1, 0.35, 0.12, lc};
xyz[] = Point{5}; // 获取点5的坐标
Point(103) = {xyz[0], xyz[1], 0.12, lc};
Line(110) = {4, 100};
Line(111) = {3, 101};
Line(112) = {6, 102};
Line(113) = {5, 103};
Line(114) = {103, 100};
Line(115) = {100, 101};
Line(116) = {101, 102};
Line(117) = {102, 103};
Curve Loop(118) = {115, -111, 3, 110};
Plane Surface(119) = {118};
Curve Loop(120) = {111, 116, -112, -7};
Plane Surface(121) = {120};
Curve Loop(122) = {112, 117, -113, -8};
Plane Surface(123) = {122};
Curve Loop(124) = {114, -110, 5, 113};
Plane Surface(125) = {124};
Curve Loop(126) = {115, 116, 117, 114};
Plane Surface(127) = {126};
Surface Loop(128) = {127, 119, 121, 123, 125, 11};
Volume(129) = {128};

// 当一个体积可以从一个面拉伸得到时, 通常使用Extrude命令更简单, 而无需手动创建所有点, 曲线和面:
Extrude {0, 0, 0.12} { Surface{my_new_surfs[1]}; }

// 以下命令允许手动为一些新点指定网格尺寸:
MeshSize {103, 105, 109, 102, 28, 24, 6, 5} = lc * 3;

// 最后, 我们将体积129和130组合到一个标签为" 1" , 名称为" The volume" 的物理组中:
Physical Volume("The volume", 1) = {129,130};

// 注意, 如果变换工具便于创建复杂几何, 有时生成" 平面" 几何(具有所有基本实体的显式表示)也是有用的.
//
// 这可以通过" 文件->导出" 并选择" Gmsh Unrolled GEO" 格式来实现, 或者在脚本中添加
//
// Save "file.geo_unrolled";
//
// 也可以通过在命令行中使用gmsh t2.geo -0来实现.

// Gmsh GEO教程3
// 拉伸网格, ONELAB参数, 选项

// 同样, 我们从包含第一个教程开始:
Include "t1.geo";

// 与t2.geo中一样, 我们计划沿z轴进行拉伸. 但在这里, 除了拉伸几何外, 我们还想拉伸2D网格.
// 这可以使用相同的Extrude命令完成, 但要指定元素" 层" (这里是2层, 第一层有8个细分, 第二层有2个细分, 两者的高度都是h/2):
h = 0.1;
Extrude {0,0,h} { Surface{1}; Layers{ {8,2}, {0.5,1} }; }

// 拉伸也可以使用旋转而不是平移来执行, 并且生成的网格可以重组为棱柱(这里我们只使用一层, 有7个细分).
// 所有旋转都由轴方向{0,1,0}, 轴点{-0.1,0,0.1}和旋转角度-Pi/2指定:
Extrude { {0,1,0} , {-0.1,0,0.1} , -Pi/2 } {
  Surface{28}; Layers{7}; Recombine;
}

// 使用内置几何内核时, 仅支持角度小于Pi的旋转. 要进行完整的旋转, 需要至少应用3次旋转. OpenCASCADE几何内核没有此限制.

// 注意, 平移{-2*h,0,0}和旋转{1,0,0}, {0,0.15,0.25}, Pi/2也可以组合形成" 扭曲" . 这里角度被指定为ONELAB参数,
// 使用DefineConstant语法. ONELAB参数可以在图形用户界面中交互式修改, 并且可以与连接到同一ONELAB数据库的其他代码交换:
DefineConstant[ angle = {90, Min 0, Max 120, Step 1, Name "Parameters/Twisting angle"} ];

// 详细来说, DefineConstant允许你将ONELAB参数" Parameters/Twisting angle" 的值分配给变量angle. 如果ONELAB数据库中不存在该参数,
// DefineConstant将创建它并分配默认值90. 此外, 如果在调用DefineConstant之前已经定义了变量angle,
// 则DefineConstant调用将被简单跳过. 这允许构建通用的参数化模型, 其参数可以从外部冻结--这些参数不再是" 参数" .

// 此功能的一个有趣用途是与-setnumber name value命令行开关结合使用, 该开关定义一个名为name, 值为value的变量.
// 调用gmsh t3.geo -setnumber angle 30将在DefineConstant之前定义angle, 使t3.geo成为非参数化的
// (" Parameters/Twisting angle" 将不会在ONELAB数据库中创建, 也无法在图形用户界面中修改).
out[] = Extrude { {-2*h,0,0}, {1,0,0} , {0,0.15,0.25} , angle * Pi / 180 } {
  Surface{50}; Layers{10}; Recombine;
};

// 在最后这个拉伸命令中, 我们通过使用Extrude命令的返回值(一个列表)以编程方式获取体积编号.
// 这个列表包含拉伸面的" 顶部" (在out[0]中), 新创建的体积(在out[1]中)以及侧面的标签(在out[2], out[3]等中).

// 然后我们可以定义一个新的物理体积(标签101)来组合所有基本体积:
Physical Volume(101) = {1, 2, out[1]};

// 现在我们来更改一些选项......由于所有交互式选项都可以在Gmsh的脚本语言中访问, 例如我们可以使点标签可见,
// 或者直接在输入文件中重新定义一些颜色:
Geometry.PointNumbers = 1;
Geometry.Color.Points = Orange;
General.Color.Text = White;
Mesh.Color.Points = {255, 0, 0};

// 注意, 所有颜色都可以用文字或数字定义, 即Mesh.Color.Points = Red等同于Mesh.Color.Points = {255,0,0};
// 并且注意, 与用户定义的变量一样, 选项既可以用作右侧值也可以用作左侧值, 因此以下命令将把面的颜色设置为与点的颜色相同:
Geometry.Color.Surfaces = Geometry.Color.Points;

// 你可以使用" 帮助->当前选项和工作区" 菜单查看所有选项的当前值. 要将所有选项保存到文件中, 使用" 文件->导出->Gmsh选项" .
// 要将当前选项与当前文件关联, 使用" 文件->保存模型选项" . 要将当前选项保存为所有未来Gmsh会话的默认选项, 使用" 文件->保存选项为默认值" .

// Gmsh GEO教程4
// 内置函数, 表面中的孔洞, 标注, 实体颜色

// 和往常一样, 我们从定义一些变量开始:
e1 = 4.5 * cm; e2 = 6 * cm / 2; e3 = 5 * cm / 2; cm = 1e-02;
h1 = 5 * cm; h2 = 10 * cm; h3 = 5 * cm; h4 = 2 * cm; h5 = 4.5 * cm; R1 = 1 * cm; R2 = 1.5 * cm; r = 1 * cm; Lc1 = 0.01; Lc2 = 0.003;

// 我们可以使用所有常见的数学函数(注意首字母大写), 以及一些有用的函数, 如Hypot(a, b) := Sqrt(a² + b²):
ccos = (-h5*R1 + e2 * Hypot(h5, Hypot(e2, R1))) / (h5² + e2²);
ssin = Sqrt(1 - ccos²);

// 然后我们使用这些变量定义一些点和线:
Point(1) = {-e1-e2, 0 , 0, Lc1}; Point(2) = {-e1-e2, h1 , 0, Lc1};
Point(3) = {-e3-r , h1 , 0, Lc2}; Point(4) = {-e3-r , h1+r , 0, Lc2}; Point(5) = {-e3 , h1+r , 0, Lc2}; Point(6) = {-e3 , h1+h2, 0, Lc1};
Point(7) = { e3 , h1+h2, 0, Lc1}; Point(8) = { e3 , h1+r , 0, Lc2};
Point(9) = { e3+r , h1+r , 0, Lc2}; Point(10)= { e3+r , h1 , 0, Lc2};
Point(11)= { e1+e2, h1 , 0, Lc1}; Point(12)= { e1+e2, 0 , 0, Lc1}; Point(13)= { e2 , 0 , 0, Lc1};
Point(14)= { R1 / ssin, h5+R1*ccos, 0, Lc2};
Point(15)= { 0 , h5 , 0, Lc2};
Point(16)= {-R1 / ssin, h5+R1*ccos, 0, Lc2};
Point(17)= {-e2 , 0.0 , 0, Lc1};
Point(18)= {-R2 , h1+h3 , 0, Lc2}; Point(19)= {-R2 , h1+h3+h4, 0, Lc2};
Point(20)= { 0 , h1+h3+h4, 0, Lc2}; Point(21)= { R2 , h1+h3+h4, 0, Lc2};
Point(22)= { R2 , h1+h3 , 0, Lc2}; Point(23)= { 0 , h1+h3 , 0, Lc2};
Point(24)= { 0, h1+h3+h4+R2, 0, Lc2}; Point(25)= { 0, h1+h3-R2, 0, Lc2};

Line(1) = {1 , 17};
Line(2) = {17, 16};

// Gmsh提供了除直线之外的其他曲线基元: 样条曲线, B样条曲线, 圆弧, 椭圆弧等. 这里我们定义一个新的圆弧, 从点14开始, 到点16结束, 圆心是点15:
Circle(3) = {14,15,16};

// 注意, 在Gmsh中, 圆弧应始终小于Pi. OpenCASCADE几何内核没有此限制.

// 然后我们可以定义额外的线和圆, 以及一个新的表面:
Line(4) = {14, 13}; Line(5) = {13, 12}; Line(6) = {12, 11}; Line(7) = {11, 10}; Circle(8) = {8, 9, 10}; Line(9) = {8, 7}; Line(10) = {7, 6}; Line(11) = {6, 5}; Line(13) = {3, 2}; Line(14) = {2, 1}; Circle(12) = {3, 4, 5}; Line(15) = {18, 19};
Circle(16) = {21, 20, 24}; Circle(17) = {24, 20, 19}; Circle(18) = {18, 23, 25}; Circle(19) = {25, 23, 22};
Line(20) = {21,22};
Curve Loop(21) = {17, -15, 18, 19, -20, 16};
Plane Surface(22) = {21};

// 但我们仍然需要定义外部表面. 由于这个表面有一个孔洞, 它的定义现在需要两个曲线环:
Curve Loop(23) = {11, -12, 13, 14, 1, 2, -3, 4, 5, 6, 7, -8, 9, 10};
Plane Surface(24) = {23, 21};

// 一般规则是, 如果一个表面有N个孔洞, 它由N+1个曲线环定义: 第一个环定义外部边界; 其他环定义孔洞的边界.

// 最后, 我们可以通过嵌入一个包含一些字符串的后处理视图来添加一些注释:
View "comments" {
  // 在窗口坐标中添加一个文本字符串, 距离左侧10像素, 距离底部10像素, 使用StrCat函数连接字符串:
  T2(10, -10, 0){ StrCat("Created on ", Today, " with Gmsh") };
  // 在模型坐标(0, 0.11, 0)处居中添加一个文本字符串:
  T3(0, 0.11, 0, TextAttributes("Align", "Center", "Font", "Helvetica")){
    "Hole"
  };
  // 如果一个字符串以'file://'开头, 其余部分被解释为图像文件. 对于3D标注, 模型坐标中的大小可以在@符号后以'widthxheight'的形式指定
  //(如果width或height中的一个为0, 则使用自然缩放; 如果两者都为0, 则使用图像的原始像素尺寸):
  T3(0, 0.09, 0, TextAttributes("Align", "Center")){
    "file://t4_image.png@0.01x0"
  };
  // 图像的3D方向可以通过在模型空间中指定图像底部和左侧边缘的方向来指定:
  T3(-0.01, 0.09, 0, 0){ "file://t4_image.png@0.01x0,0,0,1,0,1,0" };
  // 图像也可以以" 广告牌" 模式绘制, 即始终与相机平行, 通过使用#符号:
  T3(0, 0.12, 0, TextAttributes("Align", "Center")){
    "file://t4_image.png@0.01x0#"
  };
  // 2D标注的大小直接以像素给出:
  T2(350, -7, 0){ "file://t4_image.png@20x0" };
};

// 这个后处理视图采用" 解析" 格式, 即使用与.geo文件相同的解析器进行解释. 对于包含在网格上定义的实际场值的大型后处理数据集,
// 你应该使用MSH文件格式, 它允许高效存储连续或不连续的标量, 矢量和张量场, 或任意多项式阶数.

// 视图和几何实体可以响应双击事件, 这里是向控制台打印一些消息:
View[0].DoubleClickedCommand = "Printf(' View[0] has been double-clicked!' );";
Geometry.DoubleClickedCurveCommand = "Printf(' Curve %g has been double-clicked!' , Geometry.DoubleClickedEntityTag);";

// 我们也可以更改一些实体的颜色:
Color Grey50{ Surface{ 22 }; }
Color Purple{ Surface{ 24 }; }
Color Red{ Curve{ 1:14 }; }
Color Yellow{ Curve{ 15:20 }; }

// Gmsh GEO教程5
// 网格尺寸, 宏, 循环, 体积中的孔洞

// 我们从定义一些目标网格尺寸开始:
lcar1 = .1;
lcar2 = .0005;
lcar3 = .055;

// 如果我们想全局更改这些网格尺寸(不修改上述定义), 我们可以在命令行上使用-clscale选项(或在选项文件中使用Mesh.MeshSizeFactor)给出所有网格尺寸的全局缩放因子. 例如, 使用:
//
// > gmsh t5.geo -clscale 1
//
// 这个输入文件生成一个大约有3000个节点和14,000个四面体的网格. 使用
//
// > gmsh t5.geo -clscale 0.2
//
// 网格大约有231,000个节点和1,360,000个四面体. 你可以在图形用户界面中使用" 工具->统计信息" 菜单检查网格统计信息.
//
// 有关网格尺寸的更多信息, 请参见t10.geo.

// 我们继续定义一些描述截角立方体的基本实体:
Point(1) = {0.5,0.5,0.5,lcar2}; Point(2) = {0.5,0.5,0,lcar1};
Point(3) = {0,0.5,0.5,lcar1}; Point(4) = {0,0,0.5,lcar1};
Point(5) = {0.5,0,0.5,lcar1}; Point(6) = {0.5,0,0,lcar1};
Point(7) = {0,0.5,0,lcar1}; Point(8) = {0,1,0,lcar1};
Point(9) = {1,1,0,lcar1}; Point(11) = {0,1,1,lcar1}; Point(10) = {0,0,1,lcar1}; Point(12) = {1,1,1,lcar1};
Point(13) = {1,0,1,lcar1}; Point(14) = {1,0,0,lcar1};

Line(1) = {8,9}; Line(4) = {11,8}; Line(2) = {9,12}; Line(5) = {9,14}; Line(3) = {12,11}; Line(6) = {14,13};
Line(7) = {13,12}; Line(10) = {10,4}; Line(11) = {4,5}; Line(8) = {11,10}; Line(9) = {10,13}; Line(12) = {5,6};
Line(13) = {6,2}; Line(16) = {3,7}; Line(19) = {5,1}; Line(14) = {2,1}; Line(17) = {7,2}; Line(20) = {7,8}; Line(15) = {1,3}; Line(21) = {6,14}; Line(18) = {3,4};

Curve Loop(22) = {-11,-19,-15,-18}; Plane Surface(23) = {22};
Curve Loop(24) = {16,17,14,15}; Plane Surface(25) = {24};
Curve Loop(26) = {-17,20,1,5,-21,13}; Plane Surface(27) = {26};
Plane Surface(27) ={26};
Curve Loop(28) = {-4,-1,-2,-3}; Plane Surface(29) = {28};
Curve Loop(30) = {-7,2,-5,-6}; Plane Surface(31) = {30};
Curve Loop(32) = {6,-9,10,11,12,21}; Plane Surface(33) = {32};
Curve Loop(34) = {7,3,8,9}; Plane Surface(35) = {34};
Curve Loop(36) = {-10,18,-16,-20,4,-8}; Plane Surface(37) = {36};
Curve Loop(38) = {-14,-13,-12,19}; Plane Surface(39) = {38};

// 我们现在使用用户定义的宏来在立方体中雕刻一些孔洞, 而不是使用包含的文件:
Macro CheeseHole
  // 在以下命令中, 我们使用保留变量名'newp', 它会自动选择一个新的点标签. 与'newp'类似, 特殊变量'newc', 'newcl', 'news', 'newsl'和'newv'分别选择新的曲线,
  // 曲线环, 表面, 表面环和体积标签.
  //
  // 如果Geometry.OldNewReg设置为0, 新标签被选择为每个类别(点, 曲线, 曲线环等)的当前最高标签加1. 默认情况下, 为了向后兼容, Geometry.OldNewReg设置为1,
  // 并且只使用两个类别: 一个用于点, 一个用于其他所有类别.
  p1 = newp; Point(p1) = {x, y, z, lcar3};
  p2 = newp; Point(p2) = {x+r,y, z, lcar3};
  p3 = newp; Point(p3) = {x, y+r,z, lcar3};
  p4 = newp; Point(p4) = {x, y, z+r,lcar3};
  p5 = newp; Point(p5) = {x-r,y, z, lcar3};
  p6 = newp; Point(p6) = {x, y-r,z, lcar3}; p7 = newp; Point(p7) = {x, y, z-r,lcar3};

  c1 = newc; Circle(c1) = {p2,p1,p7}; c2 = newc; Circle(c2) = {p7,p1,p5}; c3 = newc; Circle(c3) = {p5,p1,p4}; c4 = newc; Circle(c4) = {p4,p1,p2};
  c5 = newc; Circle(c5) = {p2,p1,p3}; c6 = newc; Circle(c6) = {p3,p1,p5};
  c7 = newc; Circle(c7) = {p5,p1,p6}; c8 = newc; Circle(c8) = {p6,p1,p2}; c9 = newc; Circle(c9) = {p7,p1,p3}; c10 = newc; Circle(c10) = {p3,p1,p4}; c11 = newc; Circle(c11) = {p4,p1,p6}; c12 = newc; Circle(c12) = {p6,p1,p7};

  // 我们需要非平面表面来定义球形孔洞. 这里我们使用'Surface', 它可用于边界上有3或4条曲线的表面. 使用内置内核时, 如果所有曲线都是具有相同中心的圆弧,
  // 则创建一个球形补丁; 否则使用超限插值. 使用OpenCASCADE内核时, 'Surface'可用于具有任意数量边界曲线的情况, 并将通过它们拟合一个B样条补丁.
  l1 = newcl; Curve Loop(l1) = {c5,c10,c4};
  l2 = newcl; Curve Loop(l2) = {c9,-c5,c1};
  l3 = newcl; Curve Loop(l3) = {c12,-c8,-c1};
  l4 = newcl; Curve Loop(l4) = {c8,-c4,c11}; l5 = newcl; Curve Loop(l5) = {-c10,c6,c3};
  l6 = newcl; Curve Loop(l6) = {-c11,-c3,c7}; l7 = newcl; Curve Loop(l7) = {-c2,-c7,-c12};
  l8 = newcl; Curve Loop(l8) = {-c6,-c9,c2};

  s1 = news; Surface(s1) = {l1}; s2 = news; Surface(s2) = {l2};
  s3 = news; Surface(s3) = {l3}; s4 = news; Surface(s4) = {l4};
  s5 = news; Surface(s5) = {l5};
  s6 = news; Surface(s6) = {l6};
  s7 = news; Surface(s7) = {l7}; s8 = news; Surface(s8) = {l8};

  // 然后我们将表面环标签存储在一个列表中供以后参考(我们需要这些来定义最终的体积):
  theloops[t] = newsl; Surface Loop(theloops[t]) = {s1, s2, s3, s4, s5, s6, s7, s8};
  thehole = newv;
  Volume(thehole) = theloops[t];
Return

// 我们可以使用For循环在立方体中生成五个孔洞:
x = 0; y = 0.75; z = 0; r = 0.09;
For t In {1:5}
  x += 0.166; z += 0.166;
  // 我们调用'CheeseHole'宏:
  Call CheeseHole;
  // 我们为每个孔洞定义一个物理体积:
  Physical Volume (t) = thehole;
  // 我们还在终端上打印一些变量(注意, 由于.geo文件中的所有变量在内部都被视为浮点数, 格式字符串应只包含有效的浮点格式说明符, 如%g, %f, %e等):
  Printf("Hole %g (center = {%g,%g,%g}, radius = %g) has number %g!", t, x, y, z, r, thehole);
EndFor

// 然后我们可以定义立方体外表面的表面环:
theloops[0] = newreg; Surface Loop(theloops[0]) = {23:39:2};

// 立方体的体积(不包括5个孔洞)现在由6个表面环定义: 第一个表面环定义外表面; 除第一个之外的表面环定义孔洞.
//(同样, 要引用变量数组, 其标识符后跟方括号):
Volume(186) = {theloops[]};

// 注意, 使用OpenCASCADE几何内核的实体建模, 相同的几何可以用完全不同的方式构建: 参见t16.geo.

// 最后, 我们为离散化立方体(不包括孔洞)的元素定义一个物理体积(在For循环中已经为孔洞创建了物理组):
Physical Volume (10) = 186;

// 我们可以只使模型的一部分可见, 只对这个子集进行网格划分:
//
// Hide {:}
// Recursive Show { Volume{129}; }
// Mesh.MeshOnlyVisible=1;

// 可以使用选项全局更改网格算法:
Mesh.Algorithm = 6; // 用于2D网格的Frontal-Delaunay算法

// 它们也可以为单个表面设置, 例如
MeshAlgorithm Surface {31, 35} = 1; // 表面31和35上的MeshAdapt算法

// 要生成曲线网格并对其进行优化以产生可证明有效的曲线元素(参见A. Johnen, J.-F. Remacle和C. Geuzaine. Geometric validity of curvilinear finite elements.
// Journal of Computational Physics 233, pp. 359-372, 2013; 以及T. Toulorge, C. Geuzaine, J.-F. Remacle, J. Lambrechts.
// Robust untangling of curvilinear meshes. Journal of Computational Physics 254, pp. 8-26, 2013), 你可以取消注释以下行:
//
// Mesh.ElementOrder = 2;
// Mesh.HighOrderOptimize = 2;

// Gmsh GEO教程6
// 超限网格, 删除实体

// 让我们使用第一个教程中的几何作为本教程的基础:
lc = 1e-2;
Point(1) = {0, 0, 0, lc}; Point(2) = {.1, 0, 0, lc};
Point(3) = {.1, .3, 0, lc};
Point(4) = {0, .3, 0, lc};
Line(1) = {1, 2}; Line(2) = {3, 2}; Line(3) = {3, 4};
Line(4) = {4, 1};
Curve Loop(1) = {4, 1, -2, 3};
Plane Surface(1) = {1};

// 删除表面和左侧的线, 并将该线替换为3条新线:
Delete{ Surface{1}; Curve{4}; }
p1 = newp; Point(p1) = {-0.05, 0.05, 0, lc}; p2 = newp; Point(p2) = {-0.05, 0.1, 0, lc};
l1 = newc; Line(l1) = {1, p1}; l2 = newc; Line(l2) = {p1, p2}; l3 = newc; Line(l3) = {p2, 4};

// 创建一个表面:
Curve Loop(2) = {2, -1, l1, l2, l3, -3};
Plane Surface(1) = {-2};

// " 超限曲线" 网格约束明确指定曲线上节点的位置. 例如, 以下命令强制在曲线2上均匀放置20个节点(包括两个端点上的节点):
Transfinite Curve{2} = 20;

// 让我们在曲线l1, l2和l3的组合上总共放置20个点(注意点p1和p2由曲线共享, 所以我们不会创建6 + 6 + 10 = 22个节点, 而是20个! )
Transfinite Curve{l1} = 6;
Transfinite Curve{l2} = 6;
Transfinite Curve{l3} = 10;

// 最后, 我们在曲线1(反向)和曲线3上按照几何级数放置30个节点:
Transfinite Curve{-1, 3} = 30 Using Progression 1.2;

// " 超限表面" 网格约束在表面的参数平面中使用超限插值算法, 通过结构化网格连接边界上的节点. 如果表面的边界上有超过4个角点,
// 必须手动指定超限插值的角点:
Transfinite Surface{1} = {1, 2, 3, 4};

// 要创建四边形而不是三角形, 可以使用'Recombine'命令:
Recombine Surface{1};

// 当表面的边界上只有3或4个点时, 'Transfinite Surface'约束中可以省略角点列表:
Point(7) = {0.2, 0.2, 0, 1.0};
Point(8) = {0.2, 0.1, 0, 1.0}; Point(9) = {0.25, 0.2, 0, 1.0}; Point(10) = {0.3, 0.1, 0, 1.0};
Line(10) = {8, 10};
Line(11) = {10, 9}; Line(12) = {9, 7};
Line(13) = {7, 8};
Curve Loop(14) = {10, 11, 12, 13};
Plane Surface(15) = {14};
Transfinite Curve {10, 11, 12, 13} = 10;
Transfinite Surface{15};

// 三角形的生成方式可以通过在'Transfinite Surface'命令后附加"Left", "Right"或"Alternate"来控制. 尝试例如
//
// Transfinite Surface{15} Alternate;

// 最后, 我们对网格应用椭圆平滑器以获得更规则的网格:
Mesh.Smoothing = 100;

// Gmsh GEO教程7
// 背景网格

// 通过提供背景网格(即包含目标网格尺寸的后处理视图), 可以非常精确地指定网格尺寸.

// 合并一个包含目标网格尺寸的基于列表的后处理视图:
Merge "t7_bgmesh.pos";

// 如果后处理视图是基于模型的而不是基于列表的(即它基于实际网格), 我们需要创建一个新模型来包含几何, 以便对其进行网格划分不会破坏背景网格.
// 这里不需要, 因为视图是基于列表的, 但这样做没有坏处:
NewModel;

// 合并第一个教程的几何:
Merge "t1.geo";

// 将该视图应用为当前背景网格尺寸场:
Background Mesh View[0];

// 为了仅从背景网格计算网格尺寸, 而忽略任何其他尺寸约束, 可以设置:
Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;

// 参见t10.geo获取更多信息: 背景网格实际上是一般" 网格尺寸场" 的一种特殊情况.

// Gmsh GEO教程8
// 后处理, 图像导出和动画

// 除了创建几何和网格外, GEO脚本还可用于操作后处理数据集(在Gmsh中称为" 视图" ).

// 我们首先包含t1.geo以及一些后处理视图:
Include "t1.geo";
Include "view1.pos"; Include "view1.pos";
Include "view4.pos";

// Gmsh可以读取各种格式的后处理视图. 这里view1.pos和view4.pos文件采用Gmsh" 解析" 格式,
// 由GEO脚本解析器直接解释. 当然, 解析格式应该只用于相对较小的数据集: 对于较大的数据集, 使用例如MSH文件要高效得多.

// 然后我们设置一些通用选项:
General.Trackball = 0;
General.RotationX = 0; General.RotationY = 0; General.RotationZ = 0;
General.Color.Background = White; General.Color.Foreground = Black; General.Color.Text = Black;
General.Orthographic = 0; General.Axes = 0; General.SmallAxes = 0;

// 我们还为每个后处理视图设置一些选项:
v0 = PostProcessing.NbViews-4;
v1 = v0+1; v2 = v0+2; v3 = v0+3;

View[v0].IntervalsType = 2; View[v0].OffsetZ = 0.05; View[v0].RaiseZ = 0; View[v0].Light = 1; View[v0].ShowScale = 0;
View[v0].SmoothNormals = 1;

View[v1].IntervalsType = 1;
View[v1].ColorTable = { Green, Blue };
View[v1].NbIso = 10; View[v1].ShowScale = 0;

View[v2].Name = "Test...";
View[v2].Axes = 1;
View[v2].Color.Axes = Black;
View[v2].IntervalsType = 2; View[v2].Type = 2;
View[v2].AutoPosition = 0; View[v2].PositionX = 85;
View[v2].PositionY = 50; View[v2].Width = 200;
View[v2].Height = 130;

View[v3].Visible = 0;

// 你可以通过在图形用户界面中选择" 文件->导出" 直接保存MPEG电影. 有几个预定义的动画设置, 用于循环视图中的所有时间步, 或在视图之间循环.

// 但是脚本可用于构建更复杂的动画, 通过在运行时更改选项并重新渲染图形. 然后每个帧可以保存到磁盘作为图像, 多个帧可以编码形成电影.
// 下面是这样一个自定义动画的示例.
t = 0; // 初始步骤

// 从1到3循环num
For num In {1:3}
  View[v0].TimeStep = t; // 设置时间步
  View[v1].TimeStep = t;
  View[v2].TimeStep = t;
  View[v3].TimeStep = t;
  t = (View[v0].TimeStep < View[v0].NbTimeStep-1) ? t+1 : 0; // 递增

  View[v0].RaiseZ += 0.01/View[v0].Max * t; // 提升视图v0

  If (num == 3) // 当num == 3时调整图形大小, 创建640x480的帧
    General.GraphicsWidth = General.MenuWidth + 640;
    General.GraphicsHeight = 480;
  EndIf

  frames = 50;

  // 从1到frames循环num2
  For num2 In {1:frames}
    // 逐步旋转场景
    General.RotationX += 10;
    General.RotationY = General.RotationX / 3;
    General.RotationZ += 0.1;

    // 休眠0.01秒
    Sleep 0.01;

    // 绘制场景(可以使用'DrawForceChanged'代替以强制重建顶点数组, 例如在更改元素裁剪时)
    Draw;

    If(num == 3) // 取消以下行的注释以将每个帧保存为图像文件('Print'命令保存图形窗口; 'Sprintf'函数允许动态创建文件名):
      // Print Sprintf("t8-%g.gif", num2);
      // Print Sprintf("t8-%g.ppm", num2);
      // Print Sprintf("t8-%g.jpg", num2);
    EndIf
  EndFor

  If(num == 3) // ffmpeg:
    // 这里我们可以调用系统命令生成电影. 例如, 使用
    // System "ffmpeg -i t8-%d.jpg t8.mpg"
  EndIf
EndFor

// Gmsh GEO教程9
// 插件

// 插件可以添加到Gmsh中以扩展其功能. 例如, 后处理插件可以修改视图, 或基于先前加载的视图创建新视图.
// 几个默认插件与Gmsh静态链接, 例如Isosurface, CutPlane, CutSphere, Skin, Transform或Smooth.

// 插件可以与其他选项以相同的方式控制: 要么从图形界面(右键单击视图按钮, 然后选择'Plugins'), 要么从命令文件.

// 例如, 让我们包含一个三维标量视图:
Include "view3.pos" ;

// 然后我们为'Isosurface'插件设置一些选项(从3D标量视图提取等值面), 并运行它:
Plugin(Isosurface).Value = 0.67 ; // 等值级别
Plugin(Isosurface).View = 0 ; // 源视图是View[0]
Plugin(Isosurface).Run ; // 运行插件!

// 我们还为'CutPlane'插件设置一些选项(使用平面A*x+B*y+C*z+D=0计算3D视图的截面), 然后运行它:
Plugin(CutPlane).A = 0 ;
Plugin(CutPlane).B = 0.2 ;
Plugin(CutPlane).C = 1 ;
Plugin(CutPlane).D = 0 ;
Plugin(CutPlane).View = 0 ;
Plugin(CutPlane).Run ;

// 添加一个标题(按照惯例, 对于窗口坐标, 大于99999的值表示中心. 我们也可以使用'General.GraphicsWidth / 2', 但这只在当前窗口大小下使字符串居中):
Plugin(Annotate).Text = "A nice title" ;
Plugin(Annotate).X = 1.e5;
Plugin(Annotate).Y = 50 ;
Plugin(Annotate).Font = "Times-BoldItalic" ;
Plugin(Annotate).FontSize = 28 ;
Plugin(Annotate).Align = "Center" ;
Plugin(Annotate).View = 0 ;
Plugin(Annotate).Run ;

Plugin(Annotate).Text = "(and a small subtitle)" ;
Plugin(Annotate).Y = 70 ;
Plugin(Annotate).Font = "Times-Roman" ; Plugin(Annotate).FontSize = 12 ; Plugin(Annotate).Run ;

// 我们最后设置一些选项:
View[0].Light = 1;
View[0].IntervalsType = 1;
View[0].NbIso = 6;
View[0].SmoothNormals = 1;
View[1].IntervalsType = 2;
View[2].IntervalsType = 2;

// Gmsh GEO教程10
// 网格尺寸场

// 除了在几何的点上指定目标网格尺寸(参见t1.geo)或使用背景网格(参见t7.geo)外, 你还可以使用通用的网格尺寸" 场" .

// 让我们创建一个简单的矩形几何
lc = .15;
Point(1) = {0.0,0.0,0,lc}; Point(2) = {1,0.0,0,lc};
Point(3) = {1,1,0,lc}; Point(4) = {0,1,0,lc};
Point(5) = {0.2,.5,0,lc};

Line(1) = {1,2}; Line(2) = {2,3}; Line(3) = {3,4}; Line(4) = {4,1};
Curve Loop(5) = {1,2,3,4}; Plane Surface(6) = {5};

// 假设我们希望在曲线2和点5附近获得lc/30的网格元素尺寸, 在其他地方获得lc的尺寸. 为了实现这一点, 我们可以使用两个场:
// " Distance" 和" Threshold" . 我们首先在点5和曲线2上定义一个Distance场(Field[1]).
// 这个场返回到点5和(曲线上的100个等距点)曲线2的距离.
Field[1] = Distance;
Field[1].PointsList = {5};
Field[1].CurvesList = {2};
Field[1].Sampling = 100;

// 然后我们定义一个'Threshold'场, 它使用'Distance'场1的返回值来根据计算的距离定义元素尺寸的简单变化
//
// SizeMax -
//
// SizeMin -o----------------/
//          Point DistMin DistMax
Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = lc / 30;
Field[2].SizeMax = lc;
Field[2].DistMin = 0.15;
Field[2].DistMax = 0.5;

// 假设我们想使用空间坐标的数学函数来调整网格元素尺寸. 我们可以使用MathEval场来做到这一点:
Field[3] = MathEval;
Field[3].F = "cos(4*3.14*x) * sin(4*3.14*y) / 10 + 0.101";

// 我们也可以将MathEval与来自其他场的值结合起来. 例如, 让我们围绕点1定义一个'Distance'场
Field[4] = Distance;
Field[4].PointsList = {1};

// 然后我们可以创建一个'MathEval'场, 其函数依赖于'Distance'场4的返回值, 即取决于到点1的距离
//(这里使用三次定律, 最小元素尺寸= lc / 100)
Field[5] = MathEval;
Field[5].F = Sprintf("F4^3 + %g", lc / 100);

// 我们也可以使用'Box'场在盒子内部施加元素尺寸的阶跃变化
Field[6] = Box;
Field[6].VIn = lc / 15;
Field[6].VOut = lc;
Field[6].XMin = 0.3;
Field[6].XMax = 0.6;
Field[6].YMin = 0.3;
Field[6].YMax = 0.6;
Field[6].Thickness = 0.3;

// 还有许多其他类型的场可用: 参见参考手册获取完整列表. 你也可以直接在图形用户界面中通过在" 网格" 模块中选择" 定义->尺寸场" 来创建场.

// 让我们使用所有场的最小值作为背景网格尺寸场
Field[7] = Min;
Field[7].FieldsList = {2, 3, 5, 6};
Background Field = 7;

// 为了确定网格元素的尺寸, Gmsh在本地计算以下各项的最小值
//
// 1) 模型边界框的尺寸;
// 2) 如果设置了'Mesh.MeshSizeFromPoints', 则为在几何点上指定的网格尺寸;
// 3) 如果'Mesh.MeshSizeFromCurvature'为正数, 则为基于曲率的网格尺寸(该值指定每2*pi弧度的元素数量);
// 4) 背景网格尺寸场;
// 5) 任何每个实体的网格尺寸约束.
//
// 然后该值被约束在[Mesh.MeshSizeMin, Mesh.MeshSizeMax]区间内, 并乘以Mesh.MeshSizeFactor. 此外,
// 边界网格尺寸根据Mesh.MeshSizeExtendFromBoundary的值在表面和/或体积内部进行插值(默认值见相关设置).
//
// 当元素尺寸完全由网格尺寸场指定时(如本例所示), 通常希望设置
Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;

// 这将防止由于边界上的小网格尺寸而导致过度细化.

// 最后, 虽然默认的"Frontal-Delaunay"2D网格算法(Mesh.Algorithm = 6)通常会生成最高质量的网格,
// 但"Delaunay"算法(Mesh.Algorithm = 5)将更好地处理复杂的网格尺寸场--特别是具有大元素尺寸梯度的尺寸场:
Mesh.Algorithm = 5;

// Gmsh GEO教程11
// 非结构化四边形网格

// 在教程t3.geo和t6.geo中, 我们已经看到通过使用" Recombine" 关键字, 拉伸和超限网格可以" 重组" 为四边形, 棱柱或六面体. 非结构化网格也可以以相同的方式重组.
// 让我们定义一个带有解析网格尺寸场的简单几何:
Point(1) = {-1.25, -.5, 0}; Point(2) = {1.25, -.5, 0};
Point(3) = {1.25, 1.25, 0}; Point(4) = {-1.25, 1.25, 0};

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

Curve Loop(4) = {1, 2, 3, 4}; Plane Surface(100) = {4};

Field[1] = MathEval;
Field[1].F = "0.01*(1.0+30.*(y-x*x)*(y-x*x) + (1-x)*(1-x))";
Background Field = 1;

// 要生成四边形而不是三角形, 我们只需添加
Recombine Surface{100};

// 如果我们有多个表面, 可以使用'Recombine Surface {:};'. 另一种方法是指定全局选项"Mesh.RecombineAll = 1;".

// 默认的重组算法称为" Blossom" : 它使用最小成本完美匹配算法从三角剖分生成全四边形网格. 关于该算法的更多细节可以在以下论文中找到:
// J.-F. Remacle, J. Lambrechts, B. Seny, E. Marchandise, A. Johnen and C. Geuzaine, "Blossom-Quad: a non-uniform quadrilateral mesh generator using a minimum cost perfect matching algorithm",
// International Journal for Numerical Methods in Engineering 89, pp. 1102-1119, 2012.

// 对于更好的2D(平面)四边形网格, 你可以尝试实验性的" Frontal-Delaunay for quads" 网格算法, 这是一种三角剖分算法, 能够在几乎所有地方创建直角三角形:
// J.-F. Remacle, F. Henrotte, T. Carrier-Baudouin, E. Bechet, E. Marchandise, C. Geuzaine and T. Mouton. A frontal Delaunay quad mesh generator using the L^inf norm.
// International Journal for Numerical Methods in Engineering, 94, pp. 494-512, 2013.
// 取消以下行的注释以尝试用于四边形的Frontal-Delaunay算法:
// Mesh.Algorithm = 8;

// 如果重组所有三角形会导致形状不佳的四边形, 默认的重组算法可能会在网格中留下一些三角形. 在这种情况下, 要生成全四边形网格,
// 你可以细分得到的混合网格(使用Mesh.SubdivisionAlgorithm = 1), 或者使用全四边形重组算法,
// 它将自动执行更粗的网格, 然后进行重组, 平滑和细分. 取消以下行的注释以尝试全四边形算法:
// Mesh.RecombinationAlgorithm = 2; // 或3

// 注意, 你也可以在网格划分后显式应用重组算法和/或细分步骤, 如下所示:
// Mesh 2;
// RecombineMesh;
// Mesh.SubdivisionAlgorithm = 1;
// RefineMesh;

// Gmsh GEO教程12
// 使用复合实体的跨面片网格划分

// " 复合" 网格约束允许跨表面边界生成网格, 这在例如导入的CAD模型(如STEP)具有不需要的小特征时非常有用.

// 当给出" 复合曲线" 或" 复合表面" 网格约束时, 在网格生成时, Gmsh会:
// 1. 分别对基础基本几何实体进行网格划分;
// 2. 创建一个组合所有单个网格的离散实体;
// 3. 在这个离散实体上计算离散参数化(即分段线性映射);
// 4. 使用这个离散参数化而不是构成复合实体的基础基本实体的基础几何描述对离散实体进行网格划分;
// 5. 可选地, 在原始实体上重新分类网格元素和节点.

// 上面的步骤3只有在组合单个网格得到的网格可以重新参数化时才能执行, 即形状" 足够简单" . 如果形状不适合重新参数化,
// 你应该创建几何的完整网格, 并首先对其进行重新分类以生成适合重新参数化的面片(参见t13.geo).

// 步骤1中执行的单个实体的网格划分通常应该比所需的最终网格更精细; 这可以通过" Mesh.CompoundMeshSizeFactor" 选项控制.

// 步骤5中在基础基本实体上的可选重新分类由" Mesh.CompoundClassify" 选项控制.

lc = 0.1;
Point(1) = {0, 0, 0, lc}; Point(2) = {1, 0, 0, lc};
Point(3) = {1, 1, 0.5, lc}; Point(5) = {0.3, 0.2, 0, lc}; Point(7) = {0, 0.02, 0.02, lc}; Point(8) = {1, 0.05, 0.02, lc}; Point(4) = {0, 1, 0.4, lc}; Point(6) = {0, 0.01, 0.01, lc};
Point(9) = {1, 0.32, 0.02, lc};

Line(1) = {1, 2}; Line(2) = {2, 8}; Line(3) = {8, 9};
Line(4) = {9, 3}; Line(5) = {3, 4}; Line(6) = {4, 7};
Line(7) = {7, 6}; Line(8) = {6, 1}; Spline(9) = {7, 5, 9};
Line(10) = {6, 8};

Curve Loop(11) = {5, 6, 9, 4}; Surface(1) = {11};
Curve Loop(13) = {-9, 3, 10, 7}; Surface(5) = {13};
Curve Loop(15) = {-10, 2, 1, 8}; Surface(10) = {15};

// 将曲线2, 3和4视为单个曲线进行网格划分(即跨点6和7进行网格划分)
Compound Curve{2, 3, 4};

// 对曲线6, 7和8执行同样的操作
Compound Curve{6, 7, 8};

// 将表面1, 5和10视为单个表面进行网格划分(即跨曲线9和10进行网格划分)
Compound Surface{1, 5, 10};

// Gmsh GEO教程13
//
// 无需底层CAD模型对STL文件重新划分网格
//
//----
// 我们先合并一个需要重新划分网格的STL网格.
Merge "t13_data.stl";
// 首先, 我们通过沿着尖锐的几何特征分割原始表面, 对表面进行分类(" 着色" ). 这会创建新的离散表面, 曲线和点.
DefineConstant[
    // 两个三角形之间的角度超过此值时, 边会被视为尖锐边
    angle = {40, Min 20, Max 120, Step 1, Name "Parameters/表面检测角度"},
    // 对于复杂几何体, 有些面片可能过于复杂, 细长或过大而难以参数化; 设置以下选项将强制创建可参数化的面片:
    forceParametrizablePatches = {0, Choices{0,1}, Name "Parameters/创建保证可参数化的表面"},
    // 对于开放表面, 在分类过程中包含边界边:
    includeBoundary = 1,
    // 强制曲线按给定角度分割:
    curveAngle = 180
];
ClassifySurfaces{angle * Pi/180, includeBoundary, forceParametrizablePatches, curveAngle * Pi / 180};
// 为网格中的所有离散曲线和表面创建几何形状, 以便之后可以对它们重新划分网格
CreateGeometry;
// 在批处理模式下, 上述两个步骤可以通过执行" gmsh t13.stl -reparam 40" 来完成, 这会保存包含参数化信息的" t13.msh" , 该文件之后可用于重新划分网格.
// 注意, 如果有CAD模型(例如STEP文件, 参见t20.geo)而不是STL网格, 通常最好使用该CAD模型, 而不是使用通过对网格重新参数化创建的几何. 实际上, CAD几何通常更精确, 参数化更平滑, 能带来更高效, 更高质量的网格划分. Gmsh中的离散表面重新划分网格功能经过优化, 用于处理来自成像系统等没有CAD模型的密集STL网格; 它不太适合处理CAD工具为3D打印等生成的低质量STL三角化网格(这类网格为了减小文件大小而优化, 例如包含非常细长的三角形).
// 像往常一样创建一个体积
Surface Loop(1) = Surface{:};
Volume(1) = {1};
// 我们通过尺寸场来指定单元尺寸, 只是为了演示一下 :-)
funny = DefineNumber[0, Choices{0,1}, Name "Parameters/应用有趣的网格尺寸场? "];
Field[1] = MathEval;
If(funny)
    Field[1].F = "2*Sin((x+y)/5) + 3";
Else
    Field[1].F = "4";
EndIf
Background Field = 1;

// Gmsh GEO教程14
// 同调与上同调计算
//
//
// Gmsh中的同调计算利用模型的网格找到(相对)(上)同调空间基的代表性链. 这些代表性基链作为Gmsh的物理组存储在网格中, 每个链对应一个物理组.
// 创建示例几何
m = 0.5; // 网格尺寸
h = 2; // z方向高度
Point(1) = {0, 0, 0, m};
Point(2) = {10, 0, 0, m};
Point(3) = {10, 10, 0, m};
Point(4) = {0, 10, 0, m};
Point(5) = {4, 4, 0, m};
Point(6) = {6, 4, 0, m};
Point(7) = {6, 6, 0, m};
Point(8) = {4, 6, 0, m};
Point(9) = {2, 0, 0, m};
Point(10) = {8, 0, 0, m};
Point(11) = {2, 10, 0, m};
Point(12) = {8, 10, 0, m};
Line(1) = {1, 9};
Line(2) = {9, 10};
Line(3) = {10, 2};
Line(4) = {2, 3};
Line(5) = {3, 12};
Line(6) = {12, 11};
Line(7) = {11, 4};
Line(8) = {4, 1};
Line(9) = {5, 6};
Line(10) = {6, 7};
Line(11) = {7, 8};
Line(12) = {8, 5};
Curve Loop(13) = {6, 7, 8, 1, 2, 3, 4, 5};
Curve Loop(14) = {11, 12, 9, 10};
Plane Surface(15) = {13, 14};
e() = Extrude {0, 0, h}{ Surface{15}; };
// 创建物理组, 这些物理组用于定义同调计算的域和相对同调计算的子域.
// 整个域
Physical Volume(1) = {e(1)};
// 模型的四个" 终端"
Physical Surface(70) = {e(3)};
Physical Surface(71) = {e(5)};
Physical Surface(72) = {e(7)};
Physical Surface(73) = {e(9)};
// 整个域的表面
bnd() = Abs(Boundary{ Volume{e(1)}; });
Physical Surface(80) = bnd();
// 域表面中除去四个终端的部分
bnd() -= {e(3), e(5), e(7), e(9)};
Physical Surface(75) = bnd();
// 找到域模四个终端的相对同调空间基.
Homology {{1}, {70, 71, 72, 73}};
// 找到与上述基同构的同调空间基: 模非终端域表面的同调空间, 也称为细切割.
Homology {{1}, {75}};
// 找到与上述基同构的上同调空间基: 域模四个终端的上同调空间, 也称为粗切割.
Cohomology {{1}, {70, 71, 72, 73}};
// 更多示例:
// Homology {1};
// Homology;
// Homology {{1}, {80}};
// Homology {{}, {80}};
// 更多信息参见M. Pellikka, S. Suuriniemi, L. Kettunen和C. Geuzaine的<有限元建模中的同调与上同调计算>, 发表于<SIAM Journal on Scientific Computing>35(5), 第1195-1214页, 2013年.

// Gmsh GEO教程15
//
// 嵌入的点, 线和面
//
// 在默认情况下, Gmsh中跨几何维度的网格只有当低维实体位于高维实体的边界上时(即点, 曲线或面是体积的边界的一部分时)才是共形的.
// 嵌入约束允许强制网格与其他低维实体共形.
// 我们再次从包含第一个教程开始:
Include "t1.geo";
// 我们修改网格尺寸以生成更粗的网格
lc = lc * 4;
MeshSize {1:4} = lc;
// 我们定义一个新点
Point(5) = {0.02, 0.02, 0, lc};
// 可以使用" Point In Surface" 命令强制该点被" 嵌入" 到二维网格中:
Point{5} In Surface{1};
// 同样, 可以使用" Curve in Surface" 命令强制一条曲线嵌入到二维网格中:
Point(6) = {0.02, 0.12, 0, lc};
Point(7) = {0.04, 0.18, 0, lc};
Line(5) = {6, 7};
Curve{5} In Surface{1};
// 也可以使用" Curve/Point In Volume" 命令将点和曲线嵌入到体积中:
Extrude {0, 0, 0.1}{ Surface {1}; }
p = newp;
Point(p) = {0.07, 0.15, 0.025, lc};
Point{p} In Volume {1};
l = newc;
Point(p+1) = {0.025, 0.15, 0.025, lc};
Line(l) = {7, p+1};
Curve{l} In Volume {1};
// 最后, 还可以使用" Surface In Volume" 命令将一个面嵌入到体积中:
Point(p+2) = {0.02, 0.12, 0.05, lc};
Point(p+3) = {0.04, 0.12, 0.05, lc};
Point(p+4) = {0.04, 0.18, 0.05, lc};
Point(p+5) = {0.02, 0.18, 0.05, lc};
Line(l+1) = {p+2, p+3};
Line(l+2) = {p+3, p+4};
Line(l+3) = {p+4, p+5};
Line(l+4) = {p+5, p+2};
ll = newcl;
Curve Loop(ll) = {l+1:l+4};
s = news;
Plane Surface(s) = {ll};
Surface{s} In Volume {1};
// 注意, 使用OpenCASCADE内核(参见t16.geo)时, 当对不同维度的实体应用" BooleanFragments" 命令时, 如有必要, 低维实体将自动嵌入到高维实体中.
Physical Point("Embedded point") = {p};
Physical Curve("Embdded curve") = {l};
Physical Surface("Embedded surface") = {s};
Physical Volume("Volume") = {1};

// Gmsh GEO教程16
//
// 构造实体几何, OpenCASCADE几何内核
//
//
// 从版本3开始, Gmsh允许直接使用其他几何内核, 而不是以自底向上的方式使用Gmsh的内置几何内核构建模型. 这里我们使用OpenCASCADE内核:
SetFactory("OpenCASCADE");
// 我们使用构造实体几何来构建与t5.geo中相同的模型.
// 首先创建两个立方体:
Box(1) = {0,0,0, 1,1,1};
Box(2) = {0,0,0, 0.5,0.5,0.5};
// 应用布尔差运算创建" 立方体减去八分之一" 的形状:
BooleanDifference(3) = { Volume{1}; Delete; }{ Volume{2}; Delete; };
// OpenCASCADE的布尔运算总会创建新的实体. 在参数中添加" Delete" 可以自动删除原始实体.
// 然后创建五个球体:
x = 0 ; y = 0.75 ; z = 0 ; r = 0.09 ;
For t In {1:5}
    x += 0.166 ;
    z += 0.166 ;
    Sphere(3 + t) = {x,y,z,r};
    Physical Volume(t) = {3 + t};
EndFor
// 如果我们想要五个空孔, 我们会再次使用" BooleanDifference" . 但这里我们想要五个球形内含物, 其网格应与立方体的网格共形: 因此我们使用" BooleanFragments" , 它以共形方式对所有体积进行相交处理(不会创建重复的界面):
v() = BooleanFragments{ Volume{3}; Delete; }{ Volume{3 + 1 : 3 + 5}; Delete; };
// 当布尔运算导致实体发生简单修改, 并且如果使用" Delete" 删除原始实体, Gmsh会尝试为新实体分配相同的标签. (此行为由" Geometry.OCCBooleanPreserveNumbering" 选项控制. )
// 这里上面定义的" Physical Volume" 仍然有效, 因为五个球体(体积4, 5, 6, 7和8)会被碎片操作删除, 但会以相同的标签重新创建(尽管带有新的表面).
// 立方体的标签会改变, 因此我们需要通过程序来访问它:
Physical Volume(10) = v(0);
// 使用构造实体几何创建实体非常强大, 但在一些方面可能会带来实际问题, 例如为点设置网格尺寸或识别边界时.
// 要识别点或其他边界实体, 可以利用" PointfsOf" (更通用的" Boundary" 命令的特例)和" In BoundingBox" 命令.
lcar1 = .1;
lcar2 = .0005;
lcar3 = .055;
eps = 1e-3;
// 为所有体积的所有点分配网格尺寸:
MeshSize{ PointsOf{ Volume{:}; } } = lcar1;
// 覆盖五个球体的点的此约束:
MeshSize{ PointsOf{ Volume{3 + 1 : 3 + 5}; } } = lcar3;
// 通过几何搜索选择角点:
p() = Point In BoundingBox{0.5-eps, 0.5-eps, 0.5-eps, 0.5+eps, 0.5+eps, 0.5+eps};
MeshSize{ p() } = lcar2;
// 使用OpenCASCADE几何内核创建的其他示例可在" t18.geo" , " t19.geo" 和" t20.geo" 中找到, 也可在" examples/boolean" 目录中找到.

// Gmsh GEO教程17
// 各向异性背景网格
//
// 如t7.geo中所示, 通过提供背景网格(即包含目标网格尺寸的后处理视图)可以非常精确地指定网格尺寸.
// 这里, 背景网格表示为定义在正方形上的度量张量场. 在二维中, 应使用bamg作为二维网格生成器以启用各向异性网格.
SetFactory("OpenCASCADE");
// 创建一个正方形
Rectangle(1) = {-2, -2, 0, 4, 4};
// 合并包含目标各向异性网格尺寸的后处理视图
Merge "t17_bgmesh.pos";
// 将该视图应用为当前背景网格
Background Mesh View[0];
// 使用bamg
Mesh.SmoothRatio = 3;
Mesh.AnisoMax = 1000;
Mesh.Algorithm = 7;

// Gmsh GEO教程18
//
// 周期性网格
//
//
// 周期性网格约束可以施加在表面和曲线上.
// 让我们使用OpenCASCADE几何内核构建两个几何.
SetFactory("OpenCASCADE");
// 第一个几何非常简单: 一个单位立方体, 带有非均匀网格尺寸约束(故意设置以直观验证周期性约束是否有效):
Box(1) = {0, 0, 0, 1, 1, 1};
MeshSize {:} = 0.1;
MeshSize {1} = 0.02;
// 为了强制立方体表面2(右侧)的网格与表面1(左侧)的网格匹配, 设置以下周期性约束:
Periodic Surface {2} = {1} Translate {1, 0, 0};
// 在网格生成期间, 表面2上的网格将通过复制表面1上的网格来创建. 周期性约束可以通过" 平移" , " 旋转" 或一般" 仿射" 变换来指定.
// 可以以相同方式施加多个周期性约束:
Periodic Surface {6} = {5} Translate {0, 0, 1};
Periodic Surface {4} = {3} Translate {0, 1, 0};
// 对于更复杂的情况, 手动找到对应的表面可能很繁琐, 特别是当几何通过实体建模创建时. 让我们构建一个稍微复杂一些的几何.
// 我们从一个立方体和一些球体开始:
Box(10) = {2, 0, 0, 1, 1, 1};
x = 2-0.3; y = 0; z = 0;
Sphere(11) = {x, y, z, 0.35};
Sphere(12) = {x+1, y, z, 0.35};
Sphere(13) = {x, y+1, z, 0.35};
Sphere(14) = {x, y, z+1, 0.35};
Sphere(15) = {x+1, y+1, z, 0.35};
Sphere(16) = {x, y+1, z+1, 0.35};
Sphere(17) = {x+1, y, z+1, 0.35};
Sphere(18) = {x+1, y+1, z+1, 0.35};
// 首先对所有体积进行碎片操作, 这会使球体的部分突出到立方体外部:
v() = BooleanFragments { Volume{10}; Delete; }{ Volume{11:18}; Delete; };
// 要求OpenCASCADE使用STL网格计算实体更精确的边界框:
Geometry.OCCBoundsUseStl = 1;
// 然后检索原始立方体内边界框中的所有体积, 并删除外部的所有部分:
eps = 1e-3;
vin() = Volume In BoundingBox {2-eps,-eps,-eps, 2+1+eps,1+eps,1+eps};
v() -= vin();
Recursive Delete{ Volume{v()}; }
// 我们现在设置非均匀网格尺寸约束(再次为了直观检查结果):
MeshSize { PointsOf{ Volume{vin()}; }} = 0.1;
p() = Point In BoundingBox{2-eps, -eps, -eps, 2+eps, eps, eps};
MeshSize {p()} = 0.001;
// 我们现在自动识别几何左侧和右侧的对应表面.
// 首先获取左侧的所有表面:
Sxmin() = Surface In BoundingBox{2-eps, -eps, -eps, 2+eps, 1+eps, 1+eps};
For i In {0:#Sxmin()-1}
    // 然后获取每个左侧表面的边界框
    bb() = BoundingBox Surface { Sxmin(i) };
    // 将边界框平移到右侧, 并查找其中的表面:
    Sxmax() = Surface In BoundingBox { bb(0)-eps+1, bb(1)-eps, bb(2)-eps, bb(3)+eps+1, bb(4)+eps, bb(5)+eps };
    // 对于所有匹配项, 比较相应的边界框...
    For j In {0:#Sxmax()-1}
        bb2() = BoundingBox Surface { Sxmax(j) };
        bb2(0) -= 1;
        bb2(3) -= 1;
        // ...如果匹配, 则应用周期性约束
        If(Fabs(bb2(0)-bb(0)) < eps && Fabs(bb2(1)-bb(1)) < eps &&
           Fabs(bb2(2)-bb(2)) < eps && Fabs(bb2(3)-bb(3)) < eps &&
           Fabs(bb2(4)-bb(4)) < eps && Fabs(bb2(5)-bb(5)) < eps)
            Periodic Surface {Sxmax(j)} = {Sxmin(i)} Translate {1,0,0};
        EndIf
    EndFor
EndFor

// Gmsh GEO教程19
//
// 截面体, 倒圆角, 管道, 基于曲率的网格尺寸
//
//
// OpenCASCADE几何内核支持实体建模的几个有用功能.
SetFactory("OpenCASCADE");
// 可以通过(闭合)曲线环使用" ThruSections" 命令构造体积
Circle(1) = {0,0,0, 0.5};
Circle(2) = {0.1,0.05,1, 0.1};
Circle(3) = {-0.1,-0.1,2, 0.3};
Curve Loop(1) = 1;
Curve Loop(2) = 2;
Curve Loop(3) = 3;
ThruSections(1) = {1:3};
// 使用" Ruled ThruSections" 可以强制使用直纹面:
Circle(11) = {2+0,0,0, 0.5};
Circle(12) = {2+0.1,0.05,1, 0.1};
Curve Loop(11) = 11;
Curve Loop(12) = 12;
Circle(13) = {2-0.1,-0.1,2, 0.3};
Curve Loop(13) = 13;
Ruled ThruSections(11) = {11:13};
// 复制第一个体积, 并对其所有边进行倒圆角:
v() = Translate{4, 0, 0} { Duplicata{ Volume{1}; } };
f() = Abs(Boundary{ Volume{v(0)}; });
e() = Unique(Abs(Boundary{ Surface{f()}; }));
Fillet{v(0)}{e()}{0.1}
// OpenCASCADE还允许沿平滑路径进行一般拉伸. 让我们首先定义一条样条曲线:
nturns = 1;
npts = 20;
r = 1;
h = 1 * nturns;
For i In {0 : npts - 1}
    theta = i * 2*Pi*nturns/npts;
    Point(1000 + i) = {r * Cos(theta), r * Sin(theta), i * h/npts};
EndFor
Spline(1000) = {1000 : 1000 + npts - 1};
// 线框类似于曲线环, 但为开放的:
Wire(1000) = {1000};
// 我们定义想要沿样条拉伸的形状(一个圆盘):
Disk(1000) = {1,0,0, 0.2};
Rotate {{1, 0, 0}, {0, 0, 0}, Pi/2} { Surface{1000}; }
// 我们沿样条拉伸圆盘以创建管道:
Extrude { Surface{1000}; } Using Wire {1000}
// 我们删除源表面, 并增加子边的数量以更好地显示几何:
Delete{ Surface{1000}; }
Geometry.NumSubEdges = 1000;
// 我们可以激活基于曲率的网格单元尺寸计算(这里目标是每2*Pi弧度有20个单元):
Mesh.MeshSizeFromCurvature = 20;
// 我们可以约束最小和最大单元尺寸以保持在合理范围内(更多细节参见t10.geo):
Mesh.MeshSizeMin = 0.001;
Mesh.MeshSizeMax = 0.3;

//
// Gmsh GEO教程20
//
// STEP导入与操作, 几何分区
//
//
// OpenCASCADE几何内核允许导入STEP文件并对其进行修改. 在本教程中, 我们将加载一个STEP几何并将其分区为多个切片.
SetFactory("OpenCASCADE");
// 加载一个STEP文件(使用" ShapeFromFile" 而不是" Merge" 可以直接获取导入的最高维度实体的标签):
v() = ShapeFromFile("t20_data.step");
// 如果我们在合并STEP文件之前指定
//
// Geometry.OCCTargetUnit = "M";
//
// OpenCASCADE会将单位转换为米(而不是默认的毫米).
// 获取体积的边界框:
bbox() = BoundingBox Volume{v()};
xmin = bbox(0);
ymin = bbox(1);
zmin = bbox(2);
xmax = bbox(3);
ymax = bbox(4);
zmax = bbox(5);
// 我们想要将模型切成N个切片, 要么保留体积切片, 要么只保留切割得到的表面:
DefineConstant[
    N = {5, Min 2, Max 100, Step 1, Name "Parameters/0切片数量"}
    dir = {0, Choices{0="X", 1="Y", 2="Z"}, Name "Parameters/1方向"}
    surf = {0, Choices{0, 1}, Name "Parameters/2仅保留表面? "}
];
dx = (xmax - xmin);
dy = (ymax - ymin);
dz = (zmax - zmin);
L = (dir == 0) ? dz : dx;
H = (dir == 1) ? dz : dy;
// 创建第一个切割平面:
s() = {news};
If(dir == 0)
    Rectangle(s(0)) = {xmin, ymin, zmin, L, H};
ElseIf(dir == 1)
    Rotate{ {0, 1, 0}, {xmin, ymin, zmin}, -Pi/2 } { Surface{s(0)}; }
Else
    Rotate{ {1, 0, 0}, {xmin, ymin, zmin}, Pi/2 } { Surface{s(0)}; }
EndIf
tx = (dir == 0) ? dx / N : 0;
ty = (dir == 1) ? dy / N : 0;
tz = (dir == 2) ? dz / N : 0;
Translate{tx, ty, tz} { Surface{s(0)}; }
// 创建其他切割平面:
For i In {1:N-2}
    s() += Translate{i * tx, i * ty, i * tz} { Duplicata{ Surface{s(0)}; } };
EndFor
// 将体积与所有切割平面进行碎片操作(即相交):
BooleanFragments{ Volume{v()}; Delete; }{ Surface{s()}; Delete; }
// 现在删除所有不在体积边界上的表面(及其边界实体), 即" 突出" 到体积外部的切割平面部分:
Recursive Delete { Surface{:}; }
If(surf)
    // 如果我们只想保留表面, 在切割平面周围的边界框中检索表面...
    eps = 1e-4;
    s() = {};
    For i In {1:N-1}
        xx = (dir == 0) ? xmin : xmax;
        yy = (dir == 1) ? ymin : ymax;
        zz = (dir == 2) ? zmin : zmax;
        s() += Surface In BoundingBox
            {xmin - eps + i * tx, ymin - eps + i * ty, zmin - eps + i * tz,
             xx + eps + i * tx, yy + eps + i * ty, zz + eps + i * tz};
    EndFor
    // ...并删除所有其他实体:
    dels = Surface{:};
    dels -= s();
    Delete { Volume{:}; Surface{dels()}; Curve{:}; Point{:}; }
EndIf
// 最后, 指定全局网格尺寸:
Mesh.MeshSizeMin = 3;
Mesh.MeshSizeMax = 3;
// 要对网格而不是几何进行分区, 参见t21.geo.

import os
import struct
import math
import sys
if sys.platform == "win32" :
    import msvcrt
    msvcrt.setmode(0, os.O_BINARY)
    msvcrt.setmode(1, os.O_BINARY)
while(True):
    xyz = struct.unpack("ddd", os.read(0,24))
    if math.isnan(xyz[0]):
        break
    f = 0.001 + xyz[1]*0.009
    os.write(1,struct.pack("d",f))

import struct
import sys
import math
while(True):
    xyz = struct.unpack("ddd", sys.stdin.buffer.read(24))
    if math.isnan(xyz[0]):
        break
    f = 0.001 + xyz[1]*0.009
    sys.stdout.buffer.write(struct.pack("d",f))
    sys.stdout.flush()

#include <unistd.h>
int main(int argc, char **argv) {
  double xyz[3];
  while(read(STDIN_FILENO, &xyz, 3*sizeof(double)) == 3*sizeof(double)) {
    if (xyz[0] != xyz[0]) break; //nan
    double f = 0.001 + 0.009 * xyz[1];
    write(STDOUT_FILENO, &f, sizeof(double));
  }
  return 0;
}

#include <stdio.h>
#include <io.h>
#include <fcntl.h>
int main(int argc, char **argv) {
  double xyz[3];
  setmode(fileno(stdin),O_BINARY);
  setmode(fileno(stdout),O_BINARY);
  while(read(fileno(stdin), &xyz, 3*sizeof(double)) == 3*sizeof(double)) {
    if (xyz[0] != xyz[0])
      break;
    double f = f = 0.01 + 0.09 * xyz[1];
    write(fileno(stdout), &f, sizeof(double));
  }
}