EFG5_无网格法_无网格_

clear %设置区域 Length=10; Width=10; %在区域内设置节点 ndivl=10; ndivw=ndivl; [x,conn1,conn2,numnod]=mesh2(Length,Width,ndivl,ndivw); A=[1 2 5 10]; a1=A(4); a2=A(4); ro=-100; %% % n=figure; % plot(x(1,1:length(x)),x(2,1:length(x)),'o'); % hold on; %确定节点的影响域 dmax=2.5; xspac=Length/ndivl; yspac=Width/ndivw; dm(1,1:numnod)=dmax*xspac*ones(1,numnod); dm(2,1:numnod)=dmax*yspac*ones(1,numnod); %设置积分区域 % ndivlq=4; % ndivwq=5; ndivlq=ndivl/2; ndivwq=ndivw/2; [xc,conn,numcell,numq]=mesh2(Length,Width,ndivlq,ndivwq); %设置gauss点，权重以及每个区域的jacobian quado=4; %采用四点gauss积分公式 [gauss]=gauss2(quado); numq2=numcell*quado^2; gs=zeros(4,numq2); [gs]=egauss(xc,conn,gauss,numcell); % plot(gs(1,1:length(gs)),gs(2,1:length(gs)),'*'); %对gauss点进行循环组装k矩阵及f矩阵 k1=sparse(numnod,numnod); f=zeros(numnod,1); for gg=gs gpos=gg(1:2); weight=gg(3); jac=gg(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); for i=1:L en(i)=v(i); end fbody=(pi^2/a1^2+pi^2/a2^2)*sin(gpos(2)*pi/a1)*cos(gpos(1)*pi/a2)*(ro/(8.85*10^-12)); % fbody=-2*(gpos(1)^2-2*gpos(1)+gpos(2)^2-2*gpos(2)); f(en)=f(en)+weight*jac*fbody*phi'; k1(en,en)=k1(en,en)+sparse(weight*jac*((dphix'*dphix)+(dphiy'*dphiy))); end %处理边界条件 %找出自然边界条件上的节点 ind1=0; ind2=0; ind3=0; ind4=0; for j=1:numnod if(x(1,j)==0.0) % 左边界 ind1=ind1+1; n_l(1,ind1)=x(1,j); n_l(2,ind1)=x(2,j); end if(x(1,j)==10.0) % 右边界 ind2=ind2+1; n_r(1,ind2)=x(1,j); n_r(2,ind2)=x(2,j); end if (x(2,j)==10.0) % 上边界 ind3=ind3+1; n_u(1,ind3)=x(1,j); n_u(2,ind3)=x(2,j); end if (x(2,j)==0.0) % 下边界 ind4=ind4+1; n_d(1,ind4)=x(1,j); n_d(2,ind4)=x(2,j); end end %在自然边界条件上设置gauss点 %确定左边界gauss点 ind1=0; gauss=gauss2(quado); for i=1:length(n_l)-1 ycen=(n_l(2,i+1)+n_l(2,i))/2; jac=abs(n_l(2,i+1)-n_l(2,i))/2; for j=1:quado mark(j)=ycen-gauss(1,j)*jac; ind1=ind1+1; gs_l(1,ind1)=n_l(1,i); gs_l(2,ind1)=mark(j); gs_l(3,ind1)=gauss(2,j); gs_l(4,ind1)=jac; end end % plot(gs_l(1,1:length(gs_l)),gs_l(2,1:length(gs_l)),'*'); %画出左边界上gauss点 % 右边界gauss点 ind1=0; gauss=gauss2(quado); for i=1:length(n_r)-1 ycen=(n_r(2,i+1)+n_r(2,i))/2; jac=abs(n_r(2,i+1)-n_r(2,i))/2; for j=1:quado mark(j)=ycen-gauss(1,j)*jac; ind1=ind1+1; gs_r(1,ind1)=n_r(1,i); gs_r(2,ind1)=mark(j); gs_r(3,ind1)=gauss(2,j); gs_r(4,ind1)=jac; end end % plot(gs_r(1,1:length(gs_r)),gs_r(2,1:length(gs_r)),'*'); %画出左边界上gauss点 %确定上边界上的gauss点 ind2=0; gauss=gauss2(quado); for i=1:length(n_u)-1 ycen=(n_u(1,i+1)+n_u(1,i))/2; jac=abs(n_u(1,i+1)-n_u(1,i))/2; for j=1:quado mark(j)=ycen-gauss(1,j)*jac; ind2=ind2+1; gs_u(1,ind2)=mark(j); gs_u(2,ind2)=n_u(2,i); gs_u(3,ind2)=gauss(2,j); gs_u(4,ind2)=jac; end end % plot(gs_u(1,1:length(gs_u)),gs_u(2,1:length(gs_u)),'*'); %画出下边界 %确定下边界上的gauss点 ind2=0; gauss=gauss2(quado); for i=1:length(n_d)-1 ycen=(n_d(1,i+1)+n_d(1,i))/2; jac=abs(n_d(1,i+1)-n_d(1,i))/2; for j=1:quado mark(j)=ycen-gauss(1,j)*jac; ind2=ind2+1; gs_d(1,ind2)=mark(j); gs_d(2,ind2)=n_d(2,i); gs_d(3,ind2)=gauss(2,j); gs_d(4,ind2)=jac; end end % plot(gs_d(1,1:length(gs_d)),gs_d(2,1:length(gs_d)),'*'); %画出下边界 %分别对x轴，y轴上的gauss点进行循环并组装k2矩阵 rf=1000000; %设置罚函数因子 if 0 k_l=sparse(numnod,numnod); f_l=zeros(numnod,1); for gu1=gs_l gpos=gu1(1:2); weight=gu1(3); jac=gu1(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); for i=1:L en(i)=v(i); end % fbody=sin(m*pi*gpos(2)); % f_l(en)=f_l(en)+rf*weight*jac*fbody*phi'; % k_l(en,en)=k_l(en,en)+rf*weight*jac*(phi'*phi); end k_r=zeros(numnod,numnod); f_r=zeros(numnod,1); for gu1=gs_r gpos=gu1(1:2); weight=gu1(3); jac=gu1(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); for i=1:L en(i)=v(i); end % fbody=sin(n*pi*gpos(2)); % f_r(en)=f_r(en)+rf*weight*jac*fbody*phi'; % k_r(en,en)=k_r(en,en)+rf*weight*jac*(phi'*phi); end end k_u=sparse(numnod,numnod); for gu2=gs_u gpos=gu2(1:2); weight=gu2(3); jac=gu2(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); for i=1:L en(i)=v(i); end k_u(en,en)=k_u(en,en)+sparse(rf*weight*jac*(phi'*phi)); end k_d=sparse(numnod,numnod); for gu2=gs_d gpos=gu2(1:2); weight=gu2(3); jac=gu2(4); v=domain(gpos,x,dm,numnod); L=length(v); en=zeros(1,L); [phi,dphix,dphiy]=shape(gpos,dmax,x,v,dm); for i=1:L en(i)=v(i); end k_d(en,en)=k_d(en,en)+sparse(rf*weight*jac*(phi'*phi)); end %求解方程，解出节点函数值 k=k1+k_u+k_d; u=zeros(numnod,1); u=f'/k; % %精确解 uex=zeros(1,length(x)); en=0; for pos=x en=en+1; uex(en)=sin(pi*pos(2)/a1)*cos(pi*pos(1)/a2)*ro/(8.85*10^-12); end en=0; for i=1:ndivl+1 for j=1:ndivl+1 en=en+1; % 画出在x=5 线上的数值解及精确解 （ xxpos(j,i)=x(1,en); % 由分析xxpos可以看出，在x=5线上的结果，在结果矩阵的ndivl/2+1列上 yypos(j,i)=x(2,en); % yypos任何一列都可以做为画图的坐标 % ） uu(j,i)=u(en); uuex(j,i)=uex(en); end end xii=0.0:Length/ndivl:10; yii=10:-Length/ndivl:0; Erro= max(abs(u-uex))/max(abs(uex))