参考资料
G.R.Liu Y.T.GU著 王建明 周学军译 《无网格法理论及程序设计》
数值实现
Matlab 2019a
前情回顾
石中居士：无网格法与Matlab程序设计（8）——无网格径向基插值法zhuanlan.zhihu.com地球物理局 地震波场模拟实验室 无网格组
地球物理局 基建处 数值计算科
声明：
# 系列写作内容首先符合本人的研究需要，不会优先照顾读者体验。
# 仅供学习和参考，禁止转载。无网格径向基插值法：程序实现
1.支持域及影响域
构造无网格形函数的首要问题之一涉及确定局部支持域。在无网格法中，也要用到影响域的概念。
影响域被定义为一个场节点的影响区域，影响域的中心为该场节点。而支持域则为在一个计算点
可通过下面的方法利用图(b)所示的影响域选择节点进行插值。在构造一个计算点的无网格形函数时，如该计算点位于某场节点的影响域中，则该场节点属于构造该形函数的场节点。换句话说，如果一场节点的影响域包含该计算点，则该场节点参与构造该计算点的形函数。用影响域取代支持域具有如下优点：
- 影响域可很好地适应不规则分布的节点域。
- 影响域定义在问题域中的各个节点上，不同的节点可对应不同的影响域。由于各节点的影响域尺寸可以各不相同，一些节点的影响可能大于另一些节点，从而可避免由于节点分布的失衡对构造形函数的影响。
- 因为场节点数通常比积分点数少得多，故影响域的数量较支持域数量少，这将使计算过程效率更高。
由于上述原因，我们使用影响域。
一场节点的影响域可为任意形状，而影响域的尺寸可使用类似于之前所述的方法加以确定。在一二维域上，当采用矩形影响域时，分别利用沿
其中
注意对于大型问题，选择用于构造插值/近似的节点是很耗时的，故需采用某些特定算法，如桶算法（GR Liu，2002）和树算法（如可参见GR Liu和Liu，2003）。
2. 背景网格
对于一二维域可采用矩形或三角形的背景网格。三角形背景网格能较好地适应各种复杂几何形状的问题。然而为简单起见我们采用矩形背景网格。
3. 施加本质边界条件方法
因为RPIM形函数具有Kronecker
使用罚函数法所涉及的一个主要问题是如何合理地选择惩罚系数。根据FEM的做法其惩罚系数可确定为
其中
4. 径向基函数（RBFs）中的形状参数
在RPIM法中采用径向基函数构造形函数。子程序RPIM_ Shapefunc_2D中包含了复合2次(MQ)RBF、Gaussian(EXP)RBF和薄板样条(TPS)RBF。为简单起见，在此仅讨论MQ-RBF的结果，其他RBFs的结果可类似得到。
MQ-RBF中包含两个形状参数
5. 程序说明及数据结构
无网格法的分析过程如下：
- 生成问题域的几何形体，并采用一组场节点表示该问题域。
- 生成用于数值积分的背景网格。
- 通过两层循环组装系统矩阵；外层循环针对所有背景网格单元，内层循环针对每个单元的所有积分点。
- 施加本质边界条件。
- 利用标准Gaussian消元方程求解器求解系统方程。
- 利用后处理分析该问题的最终结果（位移和应力）。
主程序：MFree_Global.m
子程序：
- Input.m
功能：从外部文件输入数据。还用于计算应力-应变矩阵D。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
Input175_55.dat
*L,H,E,v,P,
48.00000 12.00000 0.3000E+8 0.30000 1000.00000
*numnode,unuse,
175 0
*ndivx,ndivy,
24 6
*numq,numcell,
55 40
*ndivxq,ndivyq,
10 4
*nquado,alf
4 100000000.000000
*Influ.domain:ALfs
3.0
*Field_nodes:x[]
1 0.00000 6.00000
2 0.00000 4.00000
3 0.00000 2.00000
4 0.00000 0.00000
5 0.00000 -2.00000
6 0.00000 -4.00000
7 0.00000 -6.00000
8 2.00000 6.00000
9 2.00000 4.00000
10 2.00000 2.00000
11 2.00000 0.00000
12 2.00000 -2.00000
13 2.00000 -4.00000
14 2.00000 -6.00000
15 4.00000 6.00000
16 4.00000 4.00000
17 4.00000 2.00000
18 4.00000 0.00000
19 4.00000 -2.00000
20 4.00000 -4.00000
21 4.00000 -6.00000
22 6.00000 6.00000
23 6.00000 4.00000
24 6.00000 2.00000
25 6.00000 0.00000
26 6.00000 -2.00000
27 6.00000 -4.00000
28 6.00000 -6.00000
29 8.00000 6.00000
30 8.00000 4.00000
31 8.00000 2.00000
32 8.00000 0.00000
33 8.00000 -2.00000
34 8.00000 -4.00000
35 8.00000 -6.00000
36 10.00000 6.00000
37 10.00000 4.00000
38 10.00000 2.00000
39 10.00000 0.00000
40 10.00000 -2.00000
41 10.00000 -4.00000
42 10.00000 -6.00000
43 12.00000 6.00000
44 12.00000 4.00000
45 12.00000 2.00000
46 12.00000 0.00000
47 12.00000 -2.00000
48 12.00000 -4.00000
49 12.00000 -6.00000
50 14.00000 6.00000
51 14.00000 4.00000
52 14.00000 2.00000
53 14.00000 0.00000
54 14.00000 -2.00000
55 14.00000 -4.00000
56 14.00000 -6.00000
57 16.00000 6.00000
58 16.00000 4.00000
59 16.00000 2.00000
60 16.00000 0.00000
61 16.00000 -2.00000
62 16.00000 -4.00000
63 16.00000 -6.00000
64 18.00000 6.00000
65 18.00000 4.00000
66 18.00000 2.00000
67 18.00000 0.00000
68 18.00000 -2.00000
69 18.00000 -4.00000
70 18.00000 -6.00000
71 20.00000 6.00000
72 20.00000 4.00000
73 20.00000 2.00000
74 20.00000 0.00000
75 20.00000 -2.00000
76 20.00000 -4.00000
77 20.00000 -6.00000
78 22.00000 6.00000
79 22.00000 4.00000
80 22.00000 2.00000
81 22.00000 0.00000
82 22.00000 -2.00000
83 22.00000 -4.00000
84 22.00000 -6.00000
85 24.00000 6.00000
86 24.00000 4.00000
87 24.00000 2.00000
88 24.00000 0.00000
89 24.00000 -2.00000
90 24.00000 -4.00000
91 24.00000 -6.00000
92 26.00000 6.00000
93 26.00000 4.00000
94 26.00000 2.00000
95 26.00000 0.00000
96 26.00000 -2.00000
97 26.00000 -4.00000
98 26.00000 -6.00000
99 28.00000 6.00000
100 28.00000 4.00000
101 28.00000 2.00000
102 28.00000 0.00000
103 28.00000 -2.00000
104 28.00000 -4.00000
105 28.00000 -6.00000
106 30.00000 6.00000
107 30.00000 4.00000
108 30.00000 2.00000
109 30.00000 0.00000
110 30.00000 -2.00000
111 30.00000 -4.00000
112 30.00000 -6.00000
113 32.00000 6.00000
114 32.00000 4.00000
115 32.00000 2.00000
116 32.00000 0.00000
117 32.00000 -2.00000
118 32.00000 -4.00000
119 32.00000 -6.00000
120 34.00000 6.00000
121 34.00000 4.00000
122 34.00000 2.00000
123 34.00000 0.00000
124 34.00000 -2.00000
125 34.00000 -4.00000
126 34.00000 -6.00000
127 36.00000 6.00000
128 36.00000 4.00000
129 36.00000 2.00000
130 36.00000 0.00000
131 36.00000 -2.00000
132 36.00000 -4.00000
133 36.00000 -6.00000
134 38.00000 6.00000
135 38.00000 4.00000
136 38.00000 2.00000
137 38.00000 0.00000
138 38.00000 -2.00000
139 38.00000 -4.00000
140 38.00000 -6.00000
141 40.00000 6.00000
142 40.00000 4.00000
143 40.00000 2.00000
144 40.00000 0.00000
145 40.00000 -2.00000
146 40.00000 -4.00000
147 40.00000 -6.00000
148 42.00000 6.00000
149 42.00000 4.00000
150 42.00000 2.00000
151 42.00000 0.00000
152 42.00000 -2.00000
153 42.00000 -4.00000
154 42.00000 -6.00000
155 44.00000 6.00000
156 44.00000 4.00000
157 44.00000 2.00000
158 44.00000 0.00000
159 44.00000 -2.00000
160 44.00000 -4.00000
161 44.00000 -6.00000
162 46.00000 6.00000
163 46.00000 4.00000
164 46.00000 2.00000
165 46.00000 0.00000
166 46.00000 -2.00000
167 46.00000 -4.00000
168 46.00000 -6.00000
169 48.00000 6.00000
170 48.00000 4.00000
171 48.00000 2.00000
172 48.00000 0.00000
173 48.00000 -2.00000
174 48.00000 -4.00000
175 48.00000 -6.00000
*Points_for_BK_cells:xc[]
1 0.00000 6.00000
2 0.00000 3.00000
3 0.00000 0.00000
4 0.00000 -3.00000
5 0.00000 -6.00000
6 4.80000 6.00000
7 4.80000 3.00000
8 4.80000 0.00000
9 4.80000 -3.00000
10 4.80000 -6.00000
11 9.60000 6.00000
12 9.60000 3.00000
13 9.60000 0.00000
14 9.60000 -3.00000
15 9.60000 -6.00000
16 14.40000 6.00000
17 14.40000 3.00000
18 14.40000 0.00000
19 14.40000 -3.00000
20 14.40000 -6.00000
21 19.20000 6.00000
22 19.20000 3.00000
23 19.20000 0.00000
24 19.20000 -3.00000
25 19.20000 -6.00000
26 24.00000 6.00000
27 24.00000 3.00000
28 24.00000 0.00000
29 24.00000 -3.00000
30 24.00000 -6.00000
31 28.80000 -6.00000
32 28.80000 3.00000
33 28.80000 0.00000
34 28.80000 -3.00000
35 28.80000 -6.00000
36 33.60000 -6.00000
37 33.60000 3.00000
38 33.60000 0.00000
39 33.60000 -3.00000
40 33.60000 -6.00000
41 38.40000 -6.00000
42 38.40000 3.00000
43 38.40000 0.00000
44 38.40000 -3.00000
45 38.40000 -6.00000
46 43.20000 -6.00000
47 43.20000 3.00000
48 43.20000 0.00000
49 43.20000 -3.00000
50 43.20000 -6.00000
51 48.00000 -6.00000
52 48.00000 3.00000
53 48.00000 0.00000
54 48.00000 -3.00000
55 48.00000 -6.00000
*Background_cells:noCell[]
1 1 2 7 6
2 2 3 8 7
3 3 4 9 8
4 4 5 10 9
5 6 7 12 11
6 7 8 13 12
7 8 9 14 13
8 9 10 15 14
9 11 12 17 16
10 12 13 18 17
11 13 14 19 18
12 14 15 20 19
13 16 17 22 21
14 17 18 23 22
15 18 19 24 23
16 19 20 25 24
17 21 22 27 26
18 22 23 28 27
19 23 24 29 28
20 24 25 30 29
21 26 27 32 31
22 27 28 33 32
23 28 29 34 33
24 29 30 35 34
25 31 32 37 36
26 32 33 38 37
27 33 34 39 38
28 34 35 40 39
29 36 37 42 41
30 37 38 43 42
31 38 39 44 43
32 39 40 45 44
33 41 42 47 46
34 42 43 48 47
35 43 44 49 48
36 44 45 50 49
37 46 47 52 51
38 47 48 53 52
39 48 49 54 53
40 49 50 55 54
*Essential_B.C.:numEBC
7
*Node,iUx,iUy,Ux,Uy
1 1 1 -0.00000E-25 -0.60000E-04
2 1 1 -0.70988E-05 -0.26667E-04
3 1 1 -0.56790E-05 -0.66667E-05
4 1 1 0.00000E-25 0.00000E-25
5 1 1 0.56790E-05 -0.66667E-05
6 1 1 0.70988E-05 -0.26667E-04
7 1 1 0.00000E-25 -0.60000E-04
*Concentrated_Natural_B.C.:numNBC
7
*Node,iTx,iTy,Tx,Ty
169 1 1 0.00000 0.0
170 1 1 0.00000 0.0
171 1 1 0.00000 0.0
172 1 1 0.00000 0.0
173 1 1 0.00000 0.0
174 1 1 0.00000 0.0
175 1 1 0.00000 0.0
*RBF shape parameters:nRBF ALFc Dc and q
1 1.0 2.0 1.03
Number of basis
3
*End of input- GaussCoefficient.m
功能：用于设置涉及标准Gauss积分的所有系数。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- CellGaussPoints.m
功能：用于设置一个背景单元中所有Gauss点并计算各Gauss点所对应的Jacobian值。现有程序需利用背景网格，其他形状（如三角形和圆形）的背景单元也可使用，使用者可对该子程序稍加修改，即可适用于其他形状的背景网格。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- SupportDomain.m
功能：用于确定一插值点的支持域以构造无网格形函数。开始计算之前给每个场节点的影响域赋值（主程序中）。搜索所有节点的影响域以得到插值所涉及的节点。如果插值点位于某个场节点的影响域中，则该场节点将被记录并被用于构造形函数的插值中。注意本程序采用矩形影响域。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- PointStiffnessMatrix.m
功能：利用
计算一积分点的刚度矩阵。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- EssentialBC.m
功能：用于施加本质边界条件。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- NaturalBC_concentrated.m
功能：用于施加集中自然边界条件。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- NaturalBC_distributed.m
功能：用于施加分布自然边界条件。针对
给定的分布自然边界条件，通过
计算节点力向量。总体力向量由各节点力向量组装而成。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- SolverBand.m
功能：该子程序为用于求解非对称带状矩阵的线性代数系统方程求解器(如Xu,1995)。事实上 无网格全局弱式法的刚度矩阵是对称的（参见4.2.2小节），但后面需使用非对称带状刚度矩阵。为避免列出太多可由标准程序库获得的标准程序，这里只提供该求解非对称带状矩阵的方程求解器。使用者可方便地调用其他对求解对称矩阵更有效的求解器。
注意，为便于理解该程序，现有程序未采用广泛使用的1维存储技术。采用与所得公式完全相同的形式将总体刚度矩阵存放于一2维数组中，并将该2维刚度矩阵直接输入该方程求解器子程序。在该子程序中须首先将原始矩阵转化为1D存储的带状矩阵。采用Gaussian消元标准方程求解器解方程。当理解了此过程后使用者也可采用功能更强大的求解器取代该求解器。
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- BandSolver.m
%%%%%% 参考《无网格法理论及程序设计》%%%%%%%%%%%%%
- GetDisplacement.m
功能：用于计算任意计算点（包括场节点）处的实际位移。如果仅考虑场节点，该子程序仅在EFG法中使用。正如之前所述，MLS近似不通过节点函数值，故由方程
求出的
其中
对于RPIM法不需利用该子程序计算场节点位移。因为RPIM形函数具有Kronecker
求出的
- GetStress.m
功能：利用
求解计算点处的应力分量。
为误差分析的需要，定义如下能量范数作为误差指标，这是因为应变或应力的精度较位移更重要
其中
子程序GetStress中用到了一个矩阵求逆子程序。该子程序GetInvasy使用的是Gauss-Jordan法。
2439
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
