参考资料
G.R.Liu. 2010. Mesh Free Methods: Moving Beyond the Finite Element Method. Second Edition.
G.R.Liu. and Y.T.Gu. 2005. An Introduction to Meshfree Methods and Their Programming
张雄.2003.无网格方法
数值实现
Matlab 2019a
地球物理局 地震波动力学研究所 无网格组
# 声明
# 欢迎批评指正,禁止转载
目录
https://zhuanlan.zhihu.com/p/174383117zhuanlan.zhihu.com代码:计算RPIM形函数
MLPG_main.m
% G.R.Liu. 2010. Mesh Free Methods: Moving Beyond the Finite Element Method. Second Edition.
% G.R.Liu. and Y.T.Gu. 2005. An Introduction to Meshfree Methods and Their Programming
% 地球物理局 地震波动力学研究所 无网格组
clear;
clc;
.......
.......
.......
[coefficient,weight] = Gauss_Coefficient(number_Gauss_points_partition,gauss);
.......
.......
.......
for node = 1:number_node
.......
.......
.......
x_quadrature_coordinates = Quadrature_Domain(x_size,y_size,x_node,y_node,x_global_boundary,x_quadrature_coordinates); % 局部正交域
.......
.......
.......
for iix = 1:number_partitions_x
.......
.......
.......
for jjy = 1:number_partitions_y
.......
.......
.......
gs = Domain_Gauss_Points(xcc,gauss,gs,dimension,...
number_Gauss_points,number_x_quadrature_coordinates,...
number_Gauss_points_partition,number_Gauss_points);
.......
.......
.......
for ie = 1:number_Gauss_points
.......
.......
.......
number_node_Shape = Support_Domain(numnode,nx,x_position,x,ds,number_node_support_domain,number_node_Shape);
.......
.......
.......
[weight,weight_xx,weight_yy] = Test_Function(domain_size,x_weight_center,x_position);
phi = RPIM_ShapeFunction_2D(gpos,x,nv,phi,...
nx,numnode,ndex,...
alfc,dc,q,nRBF,mbasis);
.......
.......
.......
end
.......
.......
.......
end
.......
.......
.......
end
.......
.......
.......
end
......
......
......
关于径向基点插值形函数原理,参见:
石中居士:无网格法理论与Matlab程序设计(6)——传统径向基点插值(RPIM)形函数zhuanlan.zhihu.comRPIM_ShapeFunction_2D.m
% G.R.Liu. 2010. Mesh Free Methods: Moving Beyond the Finite Element Method. Second Edition.
% G.R.Liu. and Y.T.Gu. 2005. An Introduction to Meshfree Methods and Their Programming
% 地球物理局 地震波动力学研究所 无网格组
function phi = RPIM_ShapeFunction_2D(gpos,x,nv,nx,ndex,alfc,dc,q,nRBF,mbasis)
rk = zeros(1,ndex + mbasis);
phi = zeros(10,ndex);
xv = zeros(nx,ndex);
rr = zeros(10,ndex + mbasis);
a = zeros(ndex + mbasis,ndex + mbasis);
g0 = zeros(ndex + mbasis,ndex + mbasis);
if nRBF == 1
rc = alfc*dc; % MQ
end
if nRBF == 2
q = alfc/(dc^2); % EXP
end
ep = 1e-20;
mg = ndex + mbasis;
g0(1:mg,1:mg) = 0;
for i = 1:ndex
nn = nv(i); % 支持域
xv(1,i) = x(1,nn); % 把局部支持域中场节点坐标保存到X中
xv(2,i) = x(2,nn);
end
%%%%%%%%%%% 组装G矩阵 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1:ndex
nn = nv(i);
rr = Compute_RadialBasis(x(1,nn),x(2,nn),xv,rr,ndex,rc,q,nRBF,mbasis);
g0(i,1:ndex) = rr(1,1:ndex);
if mbasis >= 0
g0(i,ndex+1) = 1;
g0(i,ndex+2) = x(1,nn);
g0(i,ndex+3) = x(2,nn);
g0(ndex+1,i) = 1;
g0(ndex+2,i) = x(1,nn);
g0(ndex+3,i) = x(2,nn);
end
end
%%%%%%%%%% 求解线性方程以获得形函数 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rr = Compute_RadialBasis(gpos(1),gpos(2),xv,rr,ndex,rc,q,nRBF,mbasis);
rk(1:mg) = rr(1,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
if kwji == 1
fprintf('Fail...');
pause;
end
phi(1,1:ndex) = rk(1:ndex);
%%%%%%%%%% 求解线性方程以获得dphidx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rk(1:mg) = rr(2,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
phi(2,1:ndex) = rk(1:ndex);
%%%%%%%%%% 求解线性方程以获得dphidy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rk(1:mg) = rr(3,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
phi(3,1:ndex) = rk(1:ndex);
%%%%%%%%%% 求解线性方程以获得dphidxx %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rk(1:mg) = rr(4,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
phi(4,1:ndex) = rk(1:ndex);
%%%%%%%%%% 求解线性方程以获得dphidxy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rk(1:mg) = rr(5,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
phi(5,1:ndex) = rk(1:ndex);
%%%%%%%%%% 求解线性方程以获得dphidyy %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
a(1:mg,1:mg) = g0(1:mg,1:mg);
rk(1:mg) = rr(6,1:mg);
[b,kwji] = GaussEqSolver_Sym(mg,a,rk,ep);
phi(6,1:ndex) = rk(1:ndex);
end
Compute_RadialBasis.m
% 参考《无网格法理论及程序设计》编制,源程序为Fortran语言
% 地球物理局 地震波场模拟实验室 无网格组
% 地球物理局 基建处 数值计算科
function OUTPUT = Compute_RadialBasis(x,y,xv,rk,ndex,R,q,nRBF,mbasis)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 加入线性多项式后计算径向基
% 输入:x,y,xv[],ndex,r,q,nRBF,mbasis
% nRBF = 1; MQ; 2 EXP 3 TSP
% 输出:rk[10,ndex+mbasis]
% 从1到10表示
% r,drx,dry,drdxx,drdxy,drdyy,
% drdxxx,drdxxy,drdxyy,drdyyy
rr2 = sqrt((x - xv(1,1:ndex)).^2 + (y - xv(2,1:ndex)).^2);
if nRBF == 1 % MQ
rk(1,1:ndex) = (rr2.^2 + R.^2).^q;
rk(2,1:ndex) = 2.*q.*(rr2 + R.^2).^(q-1).*(x - xv(1,1:ndex));
rk(3,1:ndex) = 2.*q.*(rr2 + R.^2).^(q-1).*(y - xv(2,1:ndex));
rk(4,1:ndex) = 2.*q.*(rr2 + R.^2).^(q-1) + 4.*(q-1).*q.*(x - xv(1,1:ndex)).^2.*(rr2 + R.^2).^(q-2);
rk(5,1:ndex) = 4.*(q-1).*q.*(x - xv(1,1:ndex)).*(y - xv(2,1:ndex)).*(rr2 + R.^2).^(q-2);
rk(6,1:ndex) = 2.*q.*(rr2 + R^2).^(q-1) + 4.*q.*(q-1).*(y - xv(2,1:ndex)).^2.*(rr2 + R.^2).^(q-2);
end
if nRBF == 2 % EXP
rk(1,1:ndex) = exp(-q*rr2^2);
rk(2,1:ndex) = -2*q*exp(-q*rr2)*(x-xv(1,1:ndex));
rk(3,1:ndex) = -2*q*exp(-q*rr2)*(y-xv(2,1:ndex));
rk(4,1:ndex) = -2*q*exp(-q*rr2) + 4*q^2*(x-xv(1,1:ndex)).^2*exp(-q*rr2);
rk(5,1:ndex) = 4*q^2*exp(-q*rr2)*(y-xv(2,1:ndex)).*(x-xv(1,1:ndex));
rk(6,1:ndex) = -2*q*exp(-q*rr2) + 4*q^2*(y-xv(2,1:ndex)).^2*exp(-q*rr2);
end
if nRBF == 3 % TSP
rk(1,1:ndex) = (rr2)^(0.5*q);
rk(2,1:ndex) = q*(x-xv(1,1:ndex))*(rr2)^(0.5*q-1);
rk(3,1:ndex) = q*(x-xv(2,1:ndex))*(rr2)^(0.5*q-1);
rk(4,1:ndex) = q*ri2^(0.5*q-1)+2*q*(0.5*q-1)*(x-xv(1,1:ndex)).^2*(rr2)^(0.5*q-2);
rk(5,1:ndex) = q*(0.5*q-1)*(x-xv(1,1:ndex)).*(y-xv(2,1:ndex))*(rr2)^(0.5*q-2);
rk(6,1:ndex) = q*ri2^(0.5*q-1)+2*q*(0.5*q-1)*(y-xv(2,1:ndex)).^2*(rr2)^(0.5*q-2);
end
if mbasis >= 0
rk(1,ndex + 1) = 1;
rk(1,ndex + 2) = x;
rk(1,ndex + 3) = y;
rk(2,ndex + 2) = 1;
rk(3,ndex + 3) = 1;
end
end
GaussEqSolver_Sym.m
function [b,kwji] = GaussEqSolver_Sym(n,a,b,ep) % ,kwji
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 采用高斯消元法求解对称线性方程ax=b
% 如果kwji=1,则无解;如果kwji=0,有解
% 输入:n,ma,a(ma,n),b(n),ep
% 输出:b,kwji
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% n:线性方程数
% ma:矩阵A行数
% A(ma,n):Ax=B的系数矩阵
% B(n):输入时为Ax=B的右边向量,输出时为结果向量
% ep:指定误差
% kwji:控制常数,当解唯一时,kwji = 0,否则kwji = 1。
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = zeros(1,n+1);
m(1:n) = (1:n);
for k = 1:n % 120
p = 0;
if abs(a(k:n,k:n))>=abs(p)
p = a(k:n,k:n);
io(k:n) = (k:n);
jo(k:n) = (k:n);
end
if (abs(p)-ep) < 0
kwji = 1; % 30
return
elseif (abs(p)-ep) == 0
kwji = 1;
return
elseif (abs(p)-ep) > 0 % 35
if (jo == k)
if(io == k)
p = 1./p;
in = n-1;
if(k == n)
b(k) = b(k)*p;
if(k == n)
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
else
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
else
a(k,k+1,in+1) = a(k,k+1,in+1)*p;
b(k) = b(k)*p;
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
else
t = a(io(k:n),k:n);
a(io(k:n),k:n) = a(k,k:n);
a(k,k:n) = t;
t = b(io(k:n));
b(io(k:n)) = b(k);
b(k) = t;
p = 1./p;
in = n - 1;
if(k == n)
b(k) = b(k)*p;
if(k == n)
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
else
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
else
a(k,k+1,in+1) = a(k,k+1,in+1)*p;
b(k) = b(k)*p;
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
end
else
t = a(1:n,jo(k:n));
a(1:n,jo(k:n)) = a(1:n,k);
a(1:n,k) = t;
t = b(io(k:n));
b(io(k:n)) = b(k);
b(k) = t;
p = 1./p;
in = n - 1;
if(k == n)
b(k) = b(k)*p;
if(k == n)
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
else
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
else
a(k,k+1,in+1) = a(k,k+1,in+1)*p;
b(k) = b(k)*p;
a(k+1:in+1,k+1:in+1) = a(k+1:in+1,k+1:in+1) - a(k+1:in+1,k).*a(k,k+1:in+1);
b(k+1:in+1) = b(k+1:in+1) - a(k+1:in+1,k).*b(k);
for i1 = 2:n
i = n + 1 - i1;
for j = i:in
b(i) = b(i) - a(i,j+1)*b(j+1);
end
end
i(1:n) = m(1:n);
a(1,i(1:n)) = b(1:n);
b(1:n) = a(1,1:n);
kwji = 0;
end
end
end
end
封面:[1]
参考
- ^https://timgsa.baidu.com/timg?image&quality=80&size=b9999_10000&sec=1598346094728&di=074bc7df2db9db0160d81c9596fa61a9&imgtype=0&src=http%3A%2F%2Fb3-q.mafengwo.net%2Fs6%2FM00%2F61%2F7F%2FwKgB4lLt3niAf8BbAAiI4kTHjvk58.jpeg