📄 poisson.html
字号:
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN" "http://www.w3.org/TR/REC-html40/loose.dtd"><html><head> <title>Description of Poisson</title> <meta name="keywords" content="Poisson"> <meta name="description" content="POISSON solve the 2-D Poisson equation"> <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1"> <meta name="generator" content="m2html © 2003 Guillaume Flandin"> <meta name="robots" content="index, follow"> <link type="text/css" rel="stylesheet" href="../../m2html.css"></head><body><a name="_top"></a><!-- # AFEM@matlab --><!-- menu.html 2_Solve --><h1>Poisson</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="box"><strong>POISSON solve the 2-D Poisson equation</strong></div><h2><a name="_synopsis"></a>SYNOPSIS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="box"><strong>function [u, energy] = Poisson(mesh, f, g_D, g_N) </strong></div><h2><a name="_description"></a>DESCRIPTION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="fragment"><pre class="comment"> POISSON solve the 2-D Poisson equation -\Delta u = f, in the current mesh with boundary conditions u = g_D on the Dirichelet boundary du/dn = g_N on the Neumann boundary USAGE [u] = Poisson(mesh, f, g_D, g_N) [u, energy] = Poisson(mesh, f, g_D, g_N) INPUT mesh: current mesh f: right side or data g_D: Dirichelet condition g_N: Neumann condition OUTPUT u: solution on the current mesh energy:energy of the discrete solution u REFERENCE Jochen Alberty, Carsten Carstensen, Stefan Funken, Remarks Around 50 Lines of MATLAB: Short Finite Element Implementation Numerical Algorithms, Volume 20, pages 117-137, 1999. NOTE It is optimized using vectorizing lauange to avoid for loops</pre></div><!-- crossreference --><h2><a name="_cross"></a>CROSS-REFERENCE INFORMATION <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2>This function calls:<ul style="list-style-image:url(../../matlabicon.gif)"></ul>This function is called by:<ul style="list-style-image:url(../../matlabicon.gif)"><li><a href="../../AFEM@matlab/1_Example/Lshape.html" class="code" title="function Lshape">Lshape</a> LSHAPE solves Poisson equation in a L-shaped domain with FEM with</li><li><a href="../../AFEM@matlab/1_Example/crack.html" class="code" title="function crack">crack</a> CRACK solves Poisson equation in a crack domain with AFEM.</li><li><a href="../../AFEM@matlab/1_Example/simple.html" class="code" title="function simple">simple</a> SIMPLE solves Poisson equation with Dirichlet boundaray condition in a</li></ul><!-- crossreference --><h2><a name="_source"></a>SOURCE CODE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="fragment"><pre><a name="_sub0" href="#_subfunctions" class="code">function [u, energy] = Poisson(mesh, f, g_D, g_N)</a><span class="comment">% POISSON solve the 2-D Poisson equation</span><span class="comment">% -\Delta u = f,</span><span class="comment">% in the current mesh with boundary conditions</span><span class="comment">% u = g_D on the Dirichelet boundary</span><span class="comment">% du/dn = g_N on the Neumann boundary</span><span class="comment">%</span><span class="comment">% USAGE</span><span class="comment">% [u] = Poisson(mesh, f, g_D, g_N)</span><span class="comment">% [u, energy] = Poisson(mesh, f, g_D, g_N)</span><span class="comment">%</span><span class="comment">% INPUT</span><span class="comment">% mesh: current mesh</span><span class="comment">% f: right side or data</span><span class="comment">% g_D: Dirichelet condition</span><span class="comment">% g_N: Neumann condition</span><span class="comment">%</span><span class="comment">% OUTPUT</span><span class="comment">% u: solution on the current mesh</span><span class="comment">% energy:energy of the discrete solution u</span><span class="comment">%</span><span class="comment">% REFERENCE</span><span class="comment">% Jochen Alberty, Carsten Carstensen, Stefan Funken,</span><span class="comment">% Remarks Around 50 Lines of MATLAB: Short Finite Element Implementation</span><span class="comment">% Numerical Algorithms, Volume 20, pages 117-137, 1999.</span><span class="comment">%</span><span class="comment">% NOTE</span><span class="comment">% It is optimized using vectorizing lauange to avoid for loops</span><span class="comment">%</span><span class="comment">% L. Chen & C. Zhang 10-08-2006</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Initialize the data</span><span class="comment">%--------------------------------------------------------------------------</span>N = size(mesh.node,1); NT = size(mesh.elem,1);A = sparse(N,N); b = sparse(N,1); u = sparse(N,1); <span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Compute vedge: edge as a vector and area of each element</span><span class="comment">%--------------------------------------------------------------------------</span>ve(:,:,1) = mesh.node(mesh.elem(:,3),:)-mesh.node(mesh.elem(:,2),:);ve(:,:,2) = mesh.node(mesh.elem(:,1),:)-mesh.node(mesh.elem(:,3),:);ve(:,:,3) = mesh.node(mesh.elem(:,2),:)-mesh.node(mesh.elem(:,1),:);area = 0.5*abs(-ve(:,1,3).*ve(:,2,2)+ve(:,2,3).*ve(:,1,2));<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Assemble Stiffness matrix</span><span class="comment">%--------------------------------------------------------------------------</span>Aij = zeros(NT,1); <span class="comment">% Aij = \int_T grad u_i * grad u_j</span><span class="keyword">for</span> i = 1:3 <span class="keyword">for</span> j = 1:3 Aij = (ve(:,1,i).*ve(:,1,j)+ve(:,2,i).*ve(:,2,j))./(4*area); A = A + sparse(mesh.elem(:,i),mesh.elem(:,j),Aij,N,N); <span class="keyword">end</span><span class="keyword">end</span><span class="comment">% More readable code is</span><span class="comment">% for j = 1 : NT</span><span class="comment">% A(elem(j,:),elem(j,:)) = A(elem(j,:),elem(j,:)) ...</span><span class="comment">% + localstiffness(node(elem(j,:),:));</span><span class="comment">% end</span><span class="comment">% with localstiffness is a function to compute lcoal stiffness matrix.</span><span class="comment">% To avoide loop for large NT, we use 'sparse' command here.</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Assemble Mass matrix</span><span class="comment">%--------------------------------------------------------------------------</span>M = sparse(mesh.elem(:,[1,1,1,2,2,2,3,3,3]), <span class="keyword">...</span> mesh.elem(:,[1,2,3,1,2,3,1,2,3]), <span class="keyword">...</span> area*[2,1,1,1,2,1,1,1,2]/12, N, N);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% More readable code is</span><span class="comment">% for j = 1 : NT</span><span class="comment">% M(elem(j,:),elem(j,:)) = area(j)/12*[2,1,1;</span><span class="comment">% 1,2,1;</span><span class="comment">% 1,1,2];</span><span class="comment">% end</span><span class="comment">% To avoide loop for large NT, we use 'sparse' command here.</span><span class="comment">% It is also very neat, by the way.</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Assemble right-hand-side</span><span class="comment">%--------------------------------------------------------------------------</span>b = M*feval(f,mesh.node);<span class="comment">% feval(f,mesh.node) is the nodal interpolation of f_I and M*f_I gives the</span><span class="comment">% right side. It is equivalent to the middle point rule</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Neumann boundary conditions</span><span class="comment">%--------------------------------------------------------------------------</span><span class="keyword">if</span> (~isempty(mesh.Neumann)) bdEdge = mesh.node(mesh.Neumann(:,1),:) <span class="keyword">...</span> - mesh.node(mesh.Neumann(:,2),:); d = sqrt(sum(bdEdge.^2,2)); mid = (mesh.node(mesh.Neumann(:,1),:) <span class="keyword">...</span> + mesh.node(mesh.Neumann(:,2),:))/2; b = b + sparse(mesh.Neumann, <span class="keyword">...</span><span class="comment"> </span> ones(size(mesh.Neumann,1),2),d.*g_N(mid)/2*[1,1],N,1); <span class="keyword">end</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Dirichlet boundary conditions</span><span class="comment">%--------------------------------------------------------------------------</span>bdNode = unique(mesh.Dirichlet);u(bdNode) = feval(g_D,mesh.node(bdNode,:));b = b - A * u; <span class="comment">% adjust the right hand side</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Solve the linear system A * U = B for the free nodes</span><span class="comment">%--------------------------------------------------------------------------</span>freeNode = find(mesh.type>0); freeNode = setdiff(freeNode, bdNode);<span class="keyword">if</span> size(freeNode)>0 u(freeNode) = A(freeNode, freeNode) \ b(freeNode);<span class="keyword">end</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Compute the energy of u</span><span class="comment">%--------------------------------------------------------------------------</span><span class="keyword">if</span> (nargout > 1) energy = full(0.5*u'*(A*u)- u'*M*feval(f,mesh.node));<span class="keyword">end</span> <span class="comment">% otherwise do not need to compute the discrete energy</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Graphic representation</span><span class="comment">%--------------------------------------------------------------------------</span>subplot(1,2,1)hold offtrisurf(mesh.elem, mesh.node(:,1), mesh.node(:,2), u', <span class="keyword">...</span> <span class="string">'FaceColor'</span>, <span class="string">'interp'</span>, <span class="string">'EdgeColor'</span>, <span class="string">'interp'</span>);view(40,15);title(<span class="string">'Solution to the Poisson equation'</span>, <span class="string">'FontSize'</span>, 14)<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% end of function POISSON</span><span class="comment">%--------------------------------------------------------------------------</span></pre></div><hr><address>Generated on Thu 02-Nov-2006 17:03:56 by <strong><a href="http://www.artefact.tk/software/matlab/m2html/" target="_parent">m2html</a></strong> © 2003</address></body></html>⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -