无网格法matlab程序

%基于局部弱式的无网格法---径向基插值法（悬臂梁问题）2011年8月29日完成
clear;
%*************************************************************************%
%步骤1：准备输入数据
%本步针对所研究的问题建模，内容包括:
%定义问题域的几何体
%生成场结点以表示该问题域。
%生成用于数值积分的背景网格。
%设置本质边界条件。
%确定参数，例如Gauss点数，影响域尺寸，RPIM的形状参数，惩罚系数等。
%*************************************************************************%
Lx=48;Ly=12;H=1;P=-1;I=1/12*Ly^3;
%图形显示几何体
X1=[0000];Y1=[00LyLy];Z1=[0HH0];C1=[0.10.30.50.7];
X2=[LxLxLxLx];Y2=[00LyLy];Z2=[0HH0];C2=[0.20.40.60.8];
X3=[0LxLx0];Y3=[00LyLy];Z3=[0000];C3=[0.30.50.70.9];
X4=[0LxLx0];Y4=[00LyLy];Z4=[HHHH];C4=[0.40.60.81];
X5=[00LxLx];Y5=[0000];Z5=[0HH0];C5=[0.70.50.30.1];
X6=[00LxLx];Y6=[LyLyLyLy];Z6=[0HH0];C6=[0.90.60.40.1];
fill3（X1Y1Z1C1X2Y2Z2C2X3Y3Z3C3X4Y4Z4C4X5Y5Z5C5X6Y6Z6C6）;
axis off
axis equal
x_min=0;x_max=Lx;
y_min=0;y_max=Ly;
%生成场结点以表示该问题域并图像显示。
nx=input（‘nx=‘）;ny=input（‘ny=‘）;
num_field_point=（nx+1）*（ny+1）;
hx=Lx/nx;hy=Ly/ny;
asx=input（‘asx=‘）;asy=input（‘asy=‘）;
rsx=asx*hx;rsy=asy*hy;       %rsxrsy为支持域的尺寸
dc=sqrt（rsx^2+rsy^2）;
aqx=input（‘aqx=‘）;aqy=input（‘aqy=‘）;
rqx=aqx*hx;rqy=aqy*hy;       %rqxrqy为局部积分域尺寸
rwx=rqx;rwy=rqy;             %rwxrwy为权函数域尺寸
field_point=zeros（num_field_point2）;
field_point（:1）=reshape（repmat（（0:hx:Lx）‘ny+11）num_field_point1）;
field_point（:2）=reshape（repmat（（0:hy:Ly）nx+11）num_field_point1）;
figure;hold on
fill（[0LxLx0][00LyLy][00.81]）
axis off
axis equal
for i=1:num_field_point
    plot（field_point（i1）field_point（i2）‘-ro‘‘LineWidth‘2...
                ‘MarkerEdgeColor‘‘k‘...
                ‘MarkerFaceColor‘‘g‘...
                ‘MarkerSize‘6）;
end
gauss_p=[-0.861136;-0.339981;0.339981;0.861136];
num_gauss_p=length（gauss_p）;
gauss_xi=reshape（repmat（gauss_p1num_gauss_p）num_gauss_p^21）;
gauss_eta=reshape（repmat（gauss_p‘num_gauss_p1）num_gauss_p^21）;
weightx=[0.3478540.6521450.6521450.347854];
weighty=weightx;
gauss_weight=reshape（weightx‘*weightynum_gauss_p^21）;
N0=[（1-gauss_xi）.*（1-gauss_eta）（1+gauss_xi）.*（1-gauss_eta）...
    （1+gauss_xi）.*（1+gauss_eta）（1-gauss_xi）.*（1+gauss_eta）]/4;
%注意局部弱式径向基点插值的结点刚度矩阵包含3部分：
%1.局部积分域内的部分。2.积分域位于求解域内部边界上的部分。3.本质边界上的部分。
%若局部积分域完全位于求解域内部，则第3部分为0。若选取的权函数域与局部积分域重合，则第2部分为0
MATRIX_K=zeros（2*num_field_point2*num_field_point）;  
MATRIX_F=zeros（2*num_field_point1）;
E=3*1e4;mu=0.3;
matrix_D=E/（1-mu^2）*[1mu0;mu10;00（1-mu）/2];
alphac=1;q=1.03;
for i=1:num_field_point
    %确定每个结点的局部积分域，积分域尺寸为rqxrqy
    local_integral_region=zeros（42）;
    left_x=field_point（i1）-rqx;
    right_x=field_point（i1）+rqx;
    down_y=field_point（i2）-rqy;
    up_y=field_point（i2）+rqy;
    left_flag=0;right_flag=0;
    if left_x<=x_min
       left_x=x_min;
       left_flag=1;
    end
    if right_x>=x_max
       right_x=x_max;
       right_flag=1;
    end
    if down_y<=y_min
       down_y=y_min;
    end
    if up_y>=y_max