2017-01-07 31 views
0

我在這裏找到了一個使用三角形元素的例子。然後,我繼續修改它生成網格的方式,用矩形元素替換三角形元素,但我不知道如何整合它們。這裏是我的版本是:用矩形元素在MATLAB中使用有限元方法逼近泊松方程

%3.6 femcode.m 

% [p,t,b] = squaregrid(m,n) % create grid of N=mn nodes to be listed in p 
% generate mesh of T=2(m-1)(n-1) right triangles in unit square 
m=6; n=5; % includes boundary nodes, mesh spacing 1/(m-1) and 1/(n-1) 
[x,y]=ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1)); % matlab forms x and y lists 
p=[x(:),y(:)]; % N by 2 matrix listing x,y coordinates of all N=mn nodes 
nelems=(m-1)*(n-1); 
t=zeros(nelems,3); 
for e=1:nelems 
    t(e) = e + floor(e/6); 
    t(e, 2) = e + floor(e/6) + 1; 
    t(e, 3) = e + floor(e/6) + 6; 
    t(e, 4) = e + floor(e/6) + 7; 
end 
% final t lists 4 node numbers of all rectangles in T by 4 matrix 
b=[1:m,m+1:m:m*n,2*m:m:m*n,m*n-m+2:m*n-1]; % bottom, left, right, top 
% b = numbers of all 2m+2n **boundary nodes** preparing for U(b)=0 

% [K,F] = assemble(p,t) % K and F for any mesh of rectangles: linear phi's 
N=size(p,1);T=nelems; % number of nodes, number of rectangles 
% p lists x,y coordinates of N nodes, t lists rectangles by 4 node numbers 
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense" 
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y) 

for e=1:T % integration over one rectangular element at a time 
    nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e 
    Pe=[ones(4,1),p(nodes,:),p(nodes,1).*p(nodes,2)]; % 4 by 4 matrix with rows=[1 x y xy] 
    Area=abs(det(Pe)); % area of triangle e = half of parallelogram area 
    C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes 
    % now compute 3 by 3 Ke and 3 by 1 Fe for element e 
    grad=C(2:3,:);Ke=Area*grad'*grad; % element matrix from slopes b,c in grad 
    Fe=Area/3; % integral of phi over triangle is volume of pyramid: f(x,y)=1 
    % multiply Fe by f at centroid for load f(x,y): one-point quadrature! 
    % centroid would be mean(p(nodes,:)) = average of 3 node coordinates 
    K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K 
    F(nodes)=F(nodes)+Fe; % add Fe to 3 components of load vector F 
end % all T element matrices and vectors now assembled into K and F 

% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0 
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b 
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F 
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K 
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb 

% Solving for the vector U will produce U(b)=0 at boundary nodes 
U=Kb\Fb % The FEM approximation is U_1 phi_1 + ... + U_N phi_N 

% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes 
% surf(p(:,1),p(:,2),0*p(:,1),U,'edgecolor','k','facecolor','interp'); 
% view(2),axis equal,colorbar 

我開始編輯的用於在三角形單元集成代碼的部分,但我不知道如何處理,或者如果它甚至是在以類似的方式完成矩形元素。

UPDATE

所以我一直試圖整合什麼Dohyun我有限的瞭解,並建議在這裏就是我現在得到:

m=12; n=12; % includes boundary nodes, mesh spacing 1/(m-1) and 1/(n-1) 
[x,y]=ndgrid((0:m-1)/(m-1),(0:n-1)/(n-1)); % matlab forms x and y lists 
p=[x(:),y(:)]; % N by 2 matrix listing x,y coordinates of all N=mn nodes 
nelems=(m-1)*(n-1); 
t=zeros(nelems,4); 
a=0; 
for e=1:nelems 
    t(e) = e + a; 
    t(e, 2) = e + a + 1; 
    t(e, 3) = e + a + m + 1; 
    t(e, 4) = e + a + m; 
    a = floor(e/n); 
end 
% final t lists 4 node numbers of all rectangles in T by 4 matrix 
b=[1:m,m+1:m:m*n,2*m:m:m*n,m*n-m+2:m*n-1]; % bottom, left, right, top 
% b = numbers of all 2m+2n **boundary nodes** preparing for U(b)=0 
N=size(p,1);T=nelems; % number of nodes, number of rectangles 
% p lists x,y coordinates of N nodes, t lists rectangles by 4 node numbers 
K=sparse(N,N); % zero matrix in sparse format: zeros(N) would be "dense" 
F=zeros(N,1); % load vector F to hold integrals of phi's times load f(x,y) 

for e=1:T % integration over one rectangular element at a time 
    nodes=t(e,:); % row of t = node numbers of the 3 corners of triangle e 
    Pe=[ones(4,1),p(nodes,:),p(nodes,1).*p(nodes,2)]; % 4 by 4 matrix with rows=[1 x y xy] 
    Area=abs(det(Pe(1:3,1:3))); % area of triangle e = half of parallelogram area 
    C=inv(Pe); % columns of C are coeffs in a+bx+cy to give phi=1,0,0 at nodes 
    % now compute 3 by 3 Ke and 3 by 1 Fe for element e 
    % grad=C(2:3,:); 
    % constantKe=Area*grad'*grad; % element matrix from slopes b,c in grad 
    for i=1:4 
    for j=1:4 
     syms x y 
     Kn = int(int(... 
       C(i,2)*C(j,2)+ ... 
       (C(i,2)*C(j,4)+C(i,4)*C(j,2))*y + ... 
       C(i,4)*C(j,4)*y^2 + ... 
       C(i,3)*C(j,3) + ... 
       (C(i,4)*C(j,3)+C(i,3)*C(j,4))*x + ... 
       C(i,4)*C(j,4)*x^2 ... 
       , x, Pe(1, 2), Pe(2, 2)), y, Pe(1, 3), Pe(3, 3)); 
     K(nodes(i),nodes(j)) = K(nodes(i),nodes(j)) + Kn; 
    end 
    end 

    Fe=Area/3; % integral of phi over triangle is volume of pyramid: f(x,y)=1 
    % multiply Fe by f at centroid for load f(x,y): one-point quadrature! 
    % centroid would be mean(p(nodes,:)) = average of 3 node coordinates 
    % K(nodes,nodes)=K(nodes,nodes)+Ke; % add Ke to 9 entries of global K 
    F(nodes)=F(nodes)+Fe; % add Fe to 4 components of load vector F 
end % all T element matrices and vectors now assembled into K and F 

% [Kb,Fb] = dirichlet(K,F,b) % assembled K was singular! K*ones(N,1)=0 
% Implement Dirichlet boundary conditions U(b)=0 at nodes in list b 
K(b,:)=0; K(:,b)=0; F(b)=0; % put zeros in boundary rows/columns of K and F 
K(b,b)=speye(length(b),length(b)); % put I into boundary submatrix of K 
Kb=K; Fb=F; % Stiffness matrix Kb (sparse format) and load vector Fb 

% Solving for the vector U will produce U(b)=0 at boundary nodes 
U=Kb\Fb; % The FEM approximation is U_1 phi_1 + ... + U_N phi_N 
U2=reshape(U',m,n); 
% Plot the FEM approximation U(x,y) with values U_1 to U_N at the nodes 
surf(U2) 

我想我可以用定積分,而不是數字的,但結果與原始程序不匹配。

Result of my program

Result of original program

回答

0

對於P1三角形元件的情況下,會發生梯度爲常數。因此整合僅僅是area*grad'*grad

但是,對於雙線性情況,梯度的內積是二階多項式。因此你需要使用數值積分。

因此,在簡單乘法中,需要另一個計算正交點基準值的循環。

此外,你的面積公式是錯誤的。搜索仿射映射。

我有github上的存儲庫,它實現了從1D到3D的Poisson方程和任意階多項式。如果你有興趣,來參觀https://github.com/dohyun64/fem_dohyun

+0

所以,糾正我,如果我錯了,但我認爲該地區應該是abs(det(Pe(1:3,1:3)))',因爲該區域的矩形元素應該是三角形元素的兩倍,我仍然可以使用Pe的前3行和列來計算。 – luckysori

+0

是的,這將工作。 – Dohyun