www.pudn.com > fem2d_nonlinear.rar > FEmesh.m, change:2013-05-28,size:1685b


function [node,dirichlet,neumann,elem]=FEmesh(xl,xr,yd,yu,M,N) 
%x1为区域左边界值(本题为-1) 
%xr为区域右边界值(本题为1) 
%yd为区域下边界值(本题为-1) 
%yu为区域上边界值(本题为1) 
%M为y轴方向剖分格数 
%N为x轴方向剖分格数 
%网格剖分,单元编号 
 
clc 
% M=10;N=10; 
m=M+1;n=N+1; 
% m=10;n=10; 
xl=-1;xr=1; 
yu=1;yd=-1; 
a=zeros(m*n);b=zeros(m*n); 
x=zeros(m*n,1);y=zeros(m*n,1); 
node=zeros(m*n,2); 
for k=1:n 
    for i=1:m 
        x(i+(k-1)*m)=xl+(k-1)*(xr-xl)/(n-1); 
        node(i+(k-1)*m,1)=x(i+(k-1)*m); 
        y(i+(k-1)*m)=yd+(i-1)*(yu-yd)/(m-1); 
        node(i+(k-1)*m,2)=y(i+(k-1)*m); 
    end 
end 
 
 
clear k i 
a(m,m)=2; 
a((n-1)*m+1,(n-1)*m+1)=2; 
a(1,1)=3; 
a(m*n,m*n)=3; 
for i=1:(n-2) 
    a(i*m+1,i*m+1)=4; 
    a((i+1)*m,(i+1)*m)=4; 
end 
clear i 
for j=2:(m-1) 
    a(j,j)=4; 
    a((n-1)*m+j,(n-1)*m+j)=4; 
end 
clear j 
for i=1:(n-2) 
    for j=2:(m-1) 
        a(i*m+j,i*m+j)=6; 
    end 
end 
clear i j 
a=sparse(a); 
% a 
clear i j 
for i=1:(m*n-1) 
    b(i,i+1)=-1; 
end 
clear i 
for i=1:(n-1)*m 
    b(i,i+m)=-1; 
end 
clear i 
for i=1:(n-1)*m-1 
    b(i,i+m+1)=-1; 
end 
 
for k=1:(n-1) 
    for h=1:(n-1) 
        b(m*k,m*h+1)=0; 
    end 
end 
% b 
clear k h 
b=sparse(b); 
b=a+b'+b; 
% b 
clear a 
a=sparse(b); 
 
clear b 
xy=[x,y]; 
gplot(a,xy,'r'); 
% Dirichlet 边界点 
for i=1:n 
    dirichlet(i)=(i-1)*m+1; 
    dirichlet(n+i)=i*m; 
end 
% Neumann 边界点 
for i=1:m 
    neumann(i)=i; 
    neumann(m-2+i)=(n-1)*m+i; 
end 
 
e=zeros(2,3); 
elem=[1,1,1]; 
for k=1:(n-1) 
    for i=1:(m-1) 
        e=[i+(k-1)*m,i+k*m,k*m+1+i;... 
            i+(k-1)*m,k*m+1+i,i+(k-1)*m+1]; 
        elem=[elem;e]; 
    end 
end 
elem(1,:)=[]; 
elem