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

```