NFE 文件

概述

为了求解一组偏微分方程，仅给出其弱形式是不够的，还需要根据方程的类型使用相应的算法。我们知道，有限元方法最终要求解的是离散后的线性代数方程组，因此，用户必须给出代数方程组的左端矩阵和右端向量的表达式。这些表达式就是由算法来实现的，并通过算法文件来给出，即这里将要介绍的 NFE 文件。从系统组件的角度来说，NFE 算法文件会被用于生成 SOLV 组件，供 GCN 文件后续调用 。

在 FEtch 系统中，我们已经把各种线性和非线性问题的常用算法放入了客户端的算法库中。用户可以直接使用其中已经设置好的算法。如果系统算法库中没有适合的算法文件可用，那么就需要用户依据对应的语法规范来设计专用的 NFE 文件了。

NFE 文件由以下 5 段信息组成：

• DEFI：全局信息声明段
• COEF：耦合变量声明段
• EQUATION：代数运算段
• SOLUTION：计算结果处理段
• FORT：Fortran 源代码插入段

其中，EQUATION 信息段和 SOLUTION 信息段均支持插入 Fortran 源代码。

本文件以 SOLUTION 信息段后空行跟一个 END 作为文件的结束标记。

下面各小节分别阐述 NFE 文件各个信息段的语法规则，每小节后面均附有该信息段的示例与说明。

DEFI 信息段

DEFI 信息段用于命名单元矩阵和荷载向量。

关键字 defi

关键字 defi 之后需要编写最多 5 行信息。

语法规则

defi
stif 名字
mass 名字
damp 名字
load 名字
(空行)

说明

1. stif 语句给出单元刚度矩阵的名字。该行是强制的，必须声明。
2. mass 语句给出单元质量矩阵的名字。如果当前物理场不含质量矩阵，可不写此行。
3. damp 语句给出单元阻尼矩阵的名字。如果当前物理场不含阻尼矩阵，可不写此行。
4. load 语句给出单元载荷向量的名字。该行是强制的，必须声明。

注：本系统的 DEFI 信息段废除了原 FEPG 系统中对 mdty、type 和 step 等语句的使用，不再读取与之相关的信息。

范例

defi  
stif s  
mass m  
damp c  
load f  
(空行)

注意：FEtch 系统规定，用户在 DEFI 信息段声明的矩阵或向量，包括单元刚度矩阵、质量矩阵、阻尼矩阵和载荷向量等，在生成 Fortran 代码段时都会在其名字前自动添加一个英文字母 "e"。例如，对于上面的范例，数组名称会变为 "es(:, :), em(:, :), ec(:, :), ef(:)" (分布式矩阵的情形) 或 "es(:, :), em(:), ec(:), ef(:)" (集中式矩阵的情形) 。这样做的目的只是为了避免与固定的组件代码命名冲突，并且保证这些名字符合 Fortran 的语法规则。因此，用户在后续插入自定义的 Fortran 代码段时，如果要直接使用或输出这些矩阵，切记在这些名字前面加上一个字母 "e"。

COEF 信息段

本信息段用来声明耦合变量，主要用于非线性迭代问题和多场耦合问题，通过它向 PDE/FBC 文件传递生成单元矩阵时可能要用到的场节点的数据。因此，COEF 信息段实际上是扮演了重要的数据通道的角色。

关键字 coef

语法规则

coef 场变量/向量名，…，场变量/向量名

范例

coef  u, u1

说明

1. coef 后面的参数可以是 "场变量名"或者 "场向量名"。这些名称都需要利用 var 或 vect 语句进行声明。具体说明详见声明语句。
2. 注意 coef 后面的变量名的顺序不是任意的。如果同时存在场向量和场变量，要求先写场向量名，后写场变量名。并且，变量名的顺序要与 PDE/FBC 文件中 coef 信息行的先后次序完全一致。
3. 由于 PDE/FBC 文件中的 coef 信息行是按场变量格式接收数据的，因此，同样是 coef 行，NFE 文件中的一个场向量名可能会对应 PDE/FBC 文件中的多个系数函数名。这与本场自由度的数量有关。
4. 系统允许主程序多传子程序少收，因此，NFE 文件中 coef 行的等价的场变量数可以多于 PDE/FBC 文件中 coef 行的场变量数，反之则不允许。

EQUATION 信息段

本信息段给出离散化后的线性代数方程组的左端矩阵和右端向量的表达式。

关键词 equation

关键词 equation 之后到下一个信息段（solution 段）之间的部分属于 EQUATION 信息段，由 5 种语句构成，即

• 声明语句
• 读写语句
• 场变量和场向量计算语句
• 左端矩阵计算语句
• 右端向量计算语句

以下分别叙述这五种语句。

声明语句

声明在本段信息中用户引入的各种变量名字的所属类型。用户可引用的类型有两种，即场变量和场向量。

关键字 var 和 vect

语法规则

var 场变量名，…，场变量名
vect 场向量名，…，场向量名

范例

var   a1, a2, b1, b2  
vect  u, u1

说明

1. var 是场变量声明语句的关键字，右边的名字规定了一个一维数组，其元素个数为节点总数。
2. vect 是场向量声明语句的关键字，右边的名字规定了一个二维数组 ，其行数与当前物理场的节点自由度相同，其列数为节点总数。
3. 一个场向量可以理解为由多个场变量组成的合集。当前的物理场有多少个自由度，一个场向量就等价于多少个场变量。
4. 多个变量名间用 "逗号" 或 "空格" 格开。
5. 用户使用的场变量和场向量必须提前进行声明，以用于开辟数组空间。
6. 同一名称在当前 NFE 文件中声明一次即可。如果前面使用了 COEF 信息段，建议在 EQUATION 信息段开头一次性声明会用到的所有的变量名称。

机理解析

在 NFE 文件中写入

vect u

系统生成的 Fortran 代码段会在 u 前自动添加前缀 e，并开辟二维动态数组空间

eu(kdgof,knode)

其中，kdgof 为节点自由度，knode 为节点总数。

同理，在 NFE 文件中写入

var u

系统也会在 u 前自动添加前缀 e，并开辟一维动态数组空间

eu(knode)

因此，为了避免出现由于使用未开辟的数组空间而导致的访问错误的情形，一定要对用到的场变量或场向量提前进行声明。

读写语句

通过读写语句用户可以按顺序从文件中读写向量（二维）数组和变量（一维）数组的值。

关键字 read 和 write

语法规则

read(s，文件名) 场变量/向量名，…，场变量/向量名
write(s，文件名) 场变量/向量名，…，场变量/向量名

范例

从 unod 文件依次读取场变量(或场向量) u 和 p。

read(s,unod) u, p

场向量 u 有三个分量，如果仅想将前两个分量写入 unodb，可编写

write(s,unodb) u(1), u(2)

说明

1. read 为读语句，write 为写语句。使用 s 是因为沿用了 FEPG 系统的命名习惯，代表顺序读写，依次读取或写入场变量或场向量。多个名称之间用 "逗号" 或 "空格" 格开。
2. 后面的参数 "场变量名" 或者 "场向量名" 都需要利用 var 或 vect 语句进行声明。因此，用户在使用读写语句操作变量之前，首先要检查一下声明语句的内容。
3. 对于场向量，读写次序为先读写完第一个自由度的所有节点上的值，再读写下一个自由度。
4. 读写语句的操作文件通常为 unod* 文件（这里用 "*" 字符代表场标识符 a、b、c … ）。该文件命名方法沿用了 FEPG 系统的命名习惯，以尽量兼容历史版本。当 * 场标识符为 a 时，通常将尾部的 * 位置留空，仅以 unod 作为文件名。

机理解析

为了加深理解，这里给出本语言 read 语句和对应的 Fortran 语句的具体转换关系。

在 NFE 文件中写入

vect u, u1  
var p

read(s,unod) u, u1(1), p

则在程序中会生成如下 Fortran 程序段：

allocatable ::eu(:,:),eu1(:,:),ep(:)
allocate(eu(kdgof,knode),eu1(kdgof,knode),ep(knode))

open(11,file='unod',form='unformatted',status='old')  
read(11)  ((eu(j,i),i=1,knode),j=1,kdgof) ,&  
(eu1(1,i),i=1,knode), &   
(ep(i),i=1,knode)  
close(11)

上述案例中，在通过 vect 和 var 对相关数组进行声明之后，read 语句开始依次读取数组的内容。write 语句的使用方法与此类似。可以看出，场向量的读写次序为按行读写，先读写完第一个自由度的所有节点上的值，再读写下一个自由度。因此，读和写操作要求彼此作用的变量的顺序必须保持一致。之前怎么存的，后面就怎么读取。尤其是在读取文件时，用户对存储文件的结构一定要有清楚的认识。

场变量和场向量计算语句

为了方便对已读入的变量和向量进行操作，本语言提供了变量和向量的计算语句。通过计算语句，用户可以将场变量和场向量作为一个数组整体，进行任何线性的计算操作。它的一般形式为：

语法规则

{场变量名}＝表达式 + {场变量名}*表达式 + … + {场变量名}*表达式
[场向量名]＝表达式 + [场向量名]*表达式 + … + [场变量名]*表达式
{场向量名(整数)}＝表达式 + {场向量名(整数)}*表达式 + … + {场向量名(整数)}*表达式

范例

{V} = {V1} - {V2}/2 
\计算场变量 V
[UE] = [U] - [U1] 
\场向量 UE 等于场向量 U 减去场向量 U1
[DU] = 0.0
\给场向量 DU 赋 0.0
{U(1)} = {U(2)}
\给出场向量 U 的第一个分量赋值

说明

1. 表达式要求满足 Fortran 的语法规范，计算所得值为实数。
2. 场变量名必须用配对花括号 "{" 和 "}" 括起来。式中所有场变量名都应已使用 var 语句声明过。
3. 场向量名必须用配对方括号 "[" 和 "]" 括起来。式中所有场向量名都应已使用 vect 语句声明过。
4. 场变量作为一个二维数组，用户可以直接操作场变量中的任何一个分量。场变量的某个分量的一般形式为： 场向量名（整数） 其中，圆括号中的整数表示对应场的第几个自由度。 一个场向量的分量其实就是一个场变量，因此，相关的计算语句继续按场变量的规则执行，即要用配对花括号 "{" 和 "}" 括起来。

左端矩阵计算语句

左端矩阵计算语句给出最终待求解的代数方程组的左端项，通常为单元刚度矩阵、单元质量矩阵和单元阻尼矩阵的线性组合。

关键词 matrix

语法规则

matrix＝±[单元矩阵名]*表达式 … ±[单元矩阵名]*表达式

范例

抛物方程的向后差分法

matrix = [s]*dt + [m]

波动方程的 Newmark 算法

matrix = [s] + [m]*a0 + [c]*a1

说明

1. 等号右边方程号中使用的单元矩阵名必须是之前在 DEFI 信息段定义过的。
2. 表达式要求满足 Fortran 的语法规范，计算所得值为实数。

右端向量计算语句

右端向量计算语句给出最终待求解的代数方程组的右端项。

关键词 forc

语法规则

forc＝±[单元矩阵名]*[单元载荷名 (或场向量名)]*表达式
   … …
   ±[单元载荷名 (或场向量名)]*表达式
   … …

范例

抛物方程的向后差分法

forc = [f]*dt+[m]*[u1]

说明

1. 等式左边的 "forc" 是本语言规定的关键字，右侧计算公式的表达方式与普通的向量代数运算一致。

2. 矩阵与向量的乘支持三种等价的表达方式，即有：

[矩阵名]*[向量名] = [矩阵名*向量名] = [矩阵名] [向量名]

SOLUTION 信息段

本段信息用于描述在求得代数方程组的解之后，程序还需要做哪些事情。例如，为后续的时间步和迭代步的计算作数据准备，或按指定格式输出结果以方便可视化。本信息段共支持 5 种语句，即

• SOLUTION 语句
• 声明语句
• 读写语句
• GiD 后处理格式输出语句
• 变量和向量计算语句

SOLUTION 语句是必需的，要求放在最前面。随后是其他语句，可按需求使用。下面分别阐述这 5 种语句的语言规则。

SOLUTION 语句

关键词 solution

语法规则

solution 场向量名

solution 语句用于声明作为代数方程组的解的场向量，以存放求解结果。

经 solution 语句声明过的场向量不必通过 vect 声明语句再次声明。

声明语句

声明语句，即 var 和 vect 语句，用于声明 SOLUTION 信息段用到的场变量和场向量，格式与 EQUATION 信息段的声明语句完全相同。因为场变量和场向量名称在当前 NFE 文件全局有效，如果有些名称已在 EQUATION 信息段声明过，此处不必再次声明。此外，前面 SOLUTION 语句声明过的场向量名，也无需再次声明。

读写语句

本段读写语句的语法规则及含义和 EQUATION 信息段的读写语句完全一致，在此不再额外阐述。

GiD 后处理格式输出语句

为了将结果按 GiD 后处理格式输出，与 FEPG 系统相比，FEtch 系统增加了一种新的写语句。

关键字 gidpost

语法规则

gidpost("标签字符串", 实数或实数变量) 场向量分量名或场变量名, … , 场向量分量名或场变量名

"标签字符串" 中对应于 GiD 软件后处理用到的场标签，它包含的字符串的数量要求与后面的场向量分量名或场变量名的数量保持一致。

"实数或实数变量" 代表输出结果的时间步编号。

用户可以通过点击客户端菜单 "后处理->res文件查看" 观察输出的 res 文件的数据格式，深入了解 gidpost 语句和 GiD 后处理文件的对应关系。

范例

稳态问题，只需输出一次结果

gidpost("T",1.0) u(1)

瞬态问题，每个时间步都要输出结果，结果与时间变量 time 相关

gidpost("ux uy",time) u(1),u(2)

说明

1. GiD 软件全面支持标量、矢量、张量等不同类型的场数据的后处理，对不同的数据类型有对应的格式要求。
2. 对于多场耦合的情形，一个节点的自由度可能是多种物理场变量的组合。只有用户自己清楚，哪些分量是真正物理意义上的标量、矢量或张量场。
3. 基于以上原因，gidpost 语句要求计算结果按照场分量的形式一个自由度一个自由度地输出，将输出的内容和顺序完全留给用户来决定。

场变量和场向量计算语句

本段场变量和场向量计算语句的语法规则及含义和 EQUATION 信息段的完全一样。除此之外，系统还提供了 %nod 和 %dof 作为对节点循环和对自由度循环的关键字。

关键字 %nod %dof

语法规则

按节点循环

%nod
循环体
%nod

按自由度循环

%dof
循环体
%dof

范例

%nod
%dof
	 [ue]=[v]-[u]
	 aa=aa+[ue]*[ue]
%dof
%nod

说明

1. 使用 %nod 或 %dof 时，符号必须成对出现，将循环体包含其中。
2. 循环变量均从 1 开始。

机理解析

为了加深理解，这里给出 %nod 和 %dof 循环体同 Fortran 语句的具体转换关系。

对于上面的范例，在程序中会生成如下 Fortran 程序段：

do i=1,knode
	do j=1,kdgof
		eue(j,i) = ev(j,i)-eu(j,i)
		aa = aa+eue(j,i)*eue(j,i)
	end do
end do

其中，knode 为节点总数，kdgof 为节点自由度。可以看出，%nod 和 %dof 分别生成了对节点和对自由度的循环体。在循环体内，生成器对场变量名执行了替换操作，替换为该场对应的数组元素，如 [v] 替换成了 ev(j,i)。这样，用户就可以方便对场变量的任一元素进行操作了，例如将其作为某个 Fortran 函数的输入或输出参数。

插入 Fortran 源代码

EQUATION 信息段和 SOLUTION 信息段均支持插入 Fortran 源代码，这与 PDE 中 DEFI 段规定的插入格式完全一致。

用户可以根据需要，在 EQUATION 信息段的 matrix 语句之前和 SOLUTION 信息段的 solution 语句之后任意位置插入 Fortran 源代码。

关键字 $c6

语法规则

插入语句的具体形式如下：

$c6 Fortran 代码

范例

简单起见，这里直接把上一个范例生成的 Fortran 代码段插入到程序中：

$c6 do i=1,knode
$c6		do j=1,kdgof
$c6			eue(j,i) = ev(j,i)-eu(j,i)
$c6			aa = aa+eue(j,i)*eue(j,i)
$c6		end do
$c6	end do

这样，我们就通过插入 Fortran 语句实现了 %nod 和 %dof 循环体的功能。该操作显然要繁琐很多，并不实用，仅供参考。

说明

注意，$c6 之后的内容会原封不动地插入到程序当中，因此用户一定要对插入语句中用到的变量和函数有充分的了解，避免出现不必要的错误。

FORT 信息段

FORT 信息段用于在指定位置嵌入一些与算法相关的 Fortran 源代码 。每段代码的前一行以符号 "@" 开始，后面给出这段代码的插入位置。编译时，生成器会自动将这段代码插入到组件代码的相应位置。

语法规则

fort
@module
Fortran 代码
@begin
Fortran 代码
@element
Fortran 代码

说明

本信息段以空行为结束行，不管插入多少段 Fortran 代码， 代码段之间都不允许有空行。

范例

对 Newmark 法的常系数变量进行定义和赋值。

fort
@begin
o=0.6
aa=0.25*(.5+o)**2
a0=1./(dt*dt*aa)
a1=o/(aa*dt)
a2=1./(aa*dt)
a3=1./(2*aa)-1.
a4=o/aa-1.
a5=dt/2*(o/aa-2.)
a6=dt*(1.-o)
a7=dt*o

调用自定义的模块。

fort  
@module  
use config

输出第一个单元的单元矩阵到指定文件。

fort
@element
! 以下ne为单元编号,k为单元数组的大小.
! es存储的是单元矩阵的转置，因此以下输出时做了转置处理.
18 format (1x,30e13.3)
if (ne .eq. 1) then
    open(21,file='matrix.txt',status='unknown')
	write(21,*) 'stif'
	do i=1,k
  		write(21,18) (es(j,i),j=1,k)
	end do 
 	write(21,*) 'load'
 	write(21,18) (ef(i),i=1,k)
 	close(21)
endif

综合案例

线性椭圆问题

线性椭圆问题是最简单的一类问题。

线性抛物问题

抛物问题是瞬态问题，它引入了时间步的更新过程，中间增加了计算结果的读写环节。

非线性椭圆问题

非线性问题是要复杂一点。它引入了迭代过程，包括临时计算结果的读写、非线性参数传递、收敛性判断等多个环节。