使用积分时在 matlab 中"Error using vertcat Dimensions of arrays being concatenated are not consistent." 2

Question

我正在编写一个关于FEM的程序，这是我的代码FEM_2D_TRI_QUA_1.m

N1=2; 
N2=2; 
N=2*N1*N2;
top=1;
bottom=0;
left=0;
right=1;
h1=(right-left)/N1;
h2=(top-bottom)/N2;
T=zeros(3,2*N1*N2);
N10=N1+1;
N20=N2+1;
for i=1:N2
for j=1:2*N1
if mod(j,2)==1
T(:,2*N1*(i-1)+j)=[(i-1)*N10+ceil(j/2);i*N10+ceil(j/2);(i-1)*N10+ceil(j/2)+1];
else
T(:,2*N1*(i-1)+j)=[(i-1)*N10+j/2+1;i*N10+j/2;i*N10+j/2+1];
end
end
end
P=zeros(2,(N1+1)*(N2+1));
for i=1:N1+1
for j=1:N2+1
P(:,(i-1)*(N1+1)+j)=[left+(i-1)*h1;bottom+(j-1)*h2];
end
end
N1_1= 2*N1;
N2_1= 2*N1;
N_1= 2*N1_1*N2_1;
h1_1=(right-left)/N1_1;
h2_1=(top-bottom)/N2_1;
N10_1=N1_1+1;
N20_1=N2_1+1;
P_b=zeros(2,(N1_1+1)*(N2_1+1));
for i=1:N1_1+1
for j=1:N2_1+1
P_b(:,(i-1)*(N1_1+1)+j)=[left+(i-1)*h1_1;bottom+(j-1)*h2_1];
end
end
boundarynodes=zeros(2,N2_1+1+N1_1+N2_1+N1_1-1);
for j=1:2*N2_1+2*N1_1
boundarynodes(1,j)=-1;
end
for i=1:N1_1+1
boundarynodes(2,i)=(i-1)*(N1_1+1)+1;
end
for i=N1_1+2:N1_1+1+N2_1
boundarynodes(2,i)=(N1_1+1-1)*(N1_1+1)+i-(N1_1);
end
for i=N1_1+1+N2_1+N1_1:-1:N1_1+N2_1+2
boundarynodes(2,i)=((-i+2*N1_1+N2_1+2)-1)*(N1_1+1)+N1_1+1;
end
for i=N1_1+N2_1+N2_1+N1_1:-1:N1_1+2+N2_1+N1_1
boundarynodes(2,i)=-i+2*N1_1+2*N2_1+2;
end
T_b=zeros(6,2*N1*N2);
c2i=@(i,j) (i-1)*(N1_1+1)+j;
for i=1:N1
for j=1:2*N2
if mod(j,2)==1
i_0=2*i-1;
j_0=2*ceil(j/2)-1;
T_b(:,2*N1*(i-1)+j)=[c2i(i_0,j_0);c2i(i_0+2,j_0);c2i(i_0,j_0+2);c2i(i_0+1,j_0);c2i(i_0+1,j_0+1);c2i(i_0,j_0+1)];
else
i_0=2*i-1;
j_0=2*(j/2+1)-1;
T_b(:,2*N1*(i-1)+j)=[c2i(i_0,j_0);c2i(i_0+2,j_0-2);c2i(i_0+2,j_0);c2i(i_0+1,j_0-1);c2i(i_0+2,j_0-1);c2i(i_0+1,j_0)];
end
end
end
A=sparse((N1_1+1)*(N2_1+1),(N1_1+1)*(N2_1+1));
for n=1:N
y_n2=P_b(2,T_b(2,n));
y_n3=P_b(2,T_b(3,n));
y_n1=P_b(2,T_b(1,n));
y_n4=P_b(2,T_b(4,n));
y_n5=P_b(2,T_b(5,n));
y_n6=P_b(2,T_b(6,n));
x_n3=P_b(1,T_b(3,n));
x_n2=P_b(1,T_b(2,n));
x_n1=P_b(1,T_b(1,n));
x_n4=P_b(1,T_b(4,n));
x_n5=P_b(1,T_b(5,n));
x_n6=P_b(1,T_b(6,n));
Y_n2=P(2,T(2,n));
Y_n3=P(2,T(3,n));
Y_n1=P(2,T(1,n));
X_n3=P(1,T(3,n));
X_n2=P(1,T(2,n));
X_n1=P(1,T(1,n));
xmin=X_n1;
xmax=xmin+h1;
J_n=(X_n2-X_n1)*(Y_n3-Y_n1)-(X_n3-X_n1)*(Y_n2-Y_n1);
x_hat=@(x,y) ((Y_n3-Y_n1)*(x-X_n1)-(X_n3-X_n1)*(y-Y_n1))/J_n;
y_hat=@(x,y) ((X_n2-X_n1)*(y-Y_n1)-(Y_n2-Y_n1)*(x-X_n1))/J_n;
if mod(n,2)==1 
ymin= Y_n1;
ymax=@(x) Y_n3+((Y_n3-Y_n2)/(X_n3-X_n2))*(x-X_n3);
else
ymin=@(x) Y_n2+((Y_n2-Y_n1)/(X_n2-X_n1))*(x-x_n2);
ymax= Y_n1;
end
for i=1:6
for j=1:6
fun=@(x,y)...
(p_i_x(x_hat(x,y),y_hat(x,y),i))*((Y_n3-Y_n1)/J_n)+...
(p_i_y(x_hat(x,y),y_hat(x,y),i))*((Y_n1-Y_n2)/J_n)*...
(p_i_x(x_hat(x,y),y_hat(x,y),j)*((Y_n3-Y_n1)/J_n)+...
(p_i_y(x_hat(x,y),y_hat(x,y),j))*((Y_n1-Y_n2)/J_n))+...
((p_i_x(x_hat(x,y),y_hat(x,y),i))*((X_n1-X_n3)/J_n)+...
(p_i_y(x_hat(x,y),y_hat(x,y),i))*((X_n2-X_n1)/J_n))*...
((p_i_x(x_hat(x,y),y_hat(x,y),j))*((X_n1-X_n3)/J_n)+...
(p_i_y(x_hat(x,y),y_hat(x,y),j))*(X_n2-X_n1)/J_n);
r=integral2(fun,xmin,xmax,ymin,ymax);
A(T_b(j,n),T_b(i,n))=A(T_b(j,n),T_b(i,n))+r;
end
end
end
A

上面的代码使用了另外两个功能代码p_i_x.m

function r=p_i_x(x,y,i)
i_x=[4*x+4*y-3;4*x-1;0;-8*x-4*y+4;4*y;-4*y];
r=i_x(i);
end

p_i_y.m

function r=p_i_y(x,y,i)
i_y=[4*x+4*y-3;0;4*y-1;-4*x;4*x;-8*y-4*x+4];
r=i_y(i);
end

我的代码在做什么？它首先将域划分为N1*N2个元素，然后我们形成关于元素和节点的信息矩阵和边界节点矩阵T，P，T_b，P_b，边界节点，然后我们想要通过嵌套for循环来组装刚度矩阵，这里我遇到了错误，这产生了

FEM_2D_TRI_QUA_1
Error using vertcat
Dimensions of arrays being concatenated are not consistent.
Error in p_i_x (line 3)
i_x=[4*x+4*y-3;4*x-1;0;-8*x-4*y+4;4*y;-4*y];
Error in
FEM_2D_TRI_QUA_1>@(x,y)(p_i_x(x_hat(x,y),y_hat(x,y),i))*((Y_n3-Y_n1)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),i))*((Y_n1-Y_n2)/J_n)*(p_i_x(x_hat(x,y),y_hat(x,y),j)*((Y_n3-Y_n1)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),j))*((Y_n1-Y_n2)/J_n))+((p_i_x(x_hat(x,y),y_hat(x,y),i))*((X_n1-X_n3)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),i))*((X_n2-X_n1)/J_n))*((p_i_x(x_hat(x,y),y_hat(x,y),j))*((X_n1-X_n3)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),j))*(X_n2-X_n1)/J_n)
Error in integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y);  NFE = NFE + 1;
Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);
Error in integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
Error in integral2 (line 106)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
Error in FEM_2D_TRI_QUA_1 (line 120)
r=integral2(fun,xmin,xmax,ymin,ymax);

根据该文件，https://ww2.mathworks.cn/help/matlab/ref/integral2.html?lang=en，函数fun必须接受两个大小相同的数组，并返回一个相应值的数组。它必须执行元素操作。但在测试了fun(0,1(后，fun正常工作，我仍然不知道哪里出了问题，有人能帮我吗，谢谢

Answer 1

我修复了它，fun@(x，y(的参数是x，y和x，y是相同大小的矩阵，所以用元素形式表示p_I_x和p_I_y中的r

function r=p_i_y(x,y,i)
if i==1
r=4*x+4*y-3;
elseif i==2
r=0*x;
elseif i==3
r=4*y-1;
elseif i==4
r=-4*x;
elseif i==5
r=4*x;
elseif i==6
r=-8*y-4*x+4;
end
end
function r=p_i_x(x,y,i)
if i==1
r=4*x+4*y-3;
elseif i==2
r=4*x-1;
elseif i==3
r=0*x;
elseif i==4
r=-8*x-4*y+4;
elseif i==5
r=4*y;
elseif i==6
r=-4*y;
end
end