%基于局部弱式的無網格法---徑向基插值法（懸臂梁問題）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
?