0

I'm coding a programme about FEM and here is my code 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

and the above code uses two other function codes 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

What is my codes doing? it first partitions the domain into N1*N2 elements and then we form the information matrix and boundarynodes matrix T,P,T_b,P_b,boundarynodes about the element and nodes, then we want to assemble the stiffness matrix via the nested for loop, here I encounter the error, which yields

 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);

according to the document,https://ww2.mathworks.cn/help/matlab/ref/integral2.html?lang=en , The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations.but after testing fun(0,1), the fun works normally, I still can't figure out where is wrong, can any one helps me, thanks you

YuerWu
  • 11
  • 1

1 Answers1

1

I fix it ,the parameters of fun @(x,y) are x,y and x,y are matrices of the same size, so express the r in p_i_x and p_i_y in elementwise form

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
YuerWu
  • 11
  • 1