📄 lshape.html

📁 C++编译指导
💻 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 Lshape</title>  <meta name="keywords" content="Lshape">  <meta name="description" content="LSHAPE solves Poisson equation in a L-shaped domain with FEM with">  <meta http-equiv="Content-Type" content="text/html; charset=iso-8859-1">  <meta name="generator" content="m2html &copy; 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 1_Example --><h1>Lshape</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="box"><strong>LSHAPE solves Poisson equation in a L-shaped domain with FEM with</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 Lshape </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"> LSHAPE solves Poisson equation in a L-shaped domain with FEM with adaptive mesh refinement.</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)"><li><a href="../../AFEM@matlab/2_Solve/Poisson.html" class="code" title="function [u, energy] = Poisson(mesh, f, g_D, g_N)">Poisson</a>	POISSON solve the 2-D Poisson equation</li><li><a href="../../AFEM@matlab/2_Solve/femerror.html" class="code" title="function [l2error, h1error] = femerror(mesh, u, u_x, u_y)">femerror</a>	</li><li><a href="../../AFEM@matlab/3_Estimate/estimate.html" class="code" title="function eta = estimate(mesh)">estimate</a>	ESTIMATE computes ZZ error estimator on each element</li><li><a href="../../AFEM@matlab/4_Refine/bisection.html" class="code" title="function [mesh,eta] = bisection(mesh,eta,theta)">bisection</a>	BISECTION	refines the triangulation using newest vertex bisection</li><li><a href="../../AFEM@matlab/4_Refine/coarsening.html" class="code" title="function [mesh,eta] = coarsening(mesh,eta,theta)">coarsening</a>	COARSENING coarsens the current mesh</li><li><a href="../../AFEM@matlab/4_Refine/uniformrefine.html" class="code" title="function [mesh] = uniformrefine(mesh)">uniformrefine</a>	UNIFORMREFINE refines the current triangulation by dividing</li><li><a href="../../AFEM@matlab/5_Mesh/getmesh.html" class="code" title="function mesh = getmesh(node,elem,Dirichlet,Neumann,type,solu)">getmesh</a>	GETMESH generates initial mesh data structure.</li></ul>This function is called by:<ul style="list-style-image:url(../../matlabicon.gif)"><li><a href="demo.html" class="code" title="function varargout = demo(varargin)">demo</a>	DEMO M-file for demo.fig</li></ul><!-- crossreference --><h2><a name="_subfunctions"></a>SUBFUNCTIONS <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><ul style="list-style-image:url(../../matlabicon.gif)"><li><a href="#_sub1" class="code">function z = f(p)</a></li><li><a href="#_sub2" class="code">function z = g_D(p)</a></li><li><a href="#_sub3" class="code">function z = g_N(p)</a></li><li><a href="#_sub4" class="code">function z = u(p)</a></li></ul><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 Lshape</a><span class="comment">% LSHAPE solves Poisson equation in a L-shaped domain with FEM with</span><span class="comment">% adaptive mesh refinement.</span><span class="comment">%</span><span class="comment">% L. Chen &amp; C. Zhang 10-20-2006</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Initialize Figure Window</span><span class="comment">%--------------------------------------------------------------------------</span>figure(1); set(gcf,<span class="string">'Units'</span>,<span class="string">'normal'</span>); set(gcf,<span class="string">'Position'</span>,[0.02,0.1,0.8,0.6]);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Parameters</span><span class="comment">%--------------------------------------------------------------------------</span>theta = 0.3; theta_c = 0.0; <span class="comment">% theta and theta_c are parameters for bisection and coarsening</span>RT = 1; CF = 8; MaxIt = 20;<span class="comment">% RT: refinement times</span><span class="comment">% CF: coarsen frequence</span><span class="comment">% MaxIt: maximum iteration number</span>N = zeros(MaxIt,1); energy = zeros(MaxIt,1); l2error = zeros(MaxIt,1);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% PreStep 0: generate initial mesh</span><span class="comment">%--------------------------------------------------------------------------</span>node = [1,0; 1,1; 0,1; -1,1; -1,0; -1,-1; 0,-1; 0,0]; <span class="comment">% nodes</span>elem = [1,2,8; 3,8,2; 8,3,5; 4,5,3; 7,8,6; 5,6,8];    <span class="comment">% elements</span>Dirichlet = [1,2; 2,3; 3,4; 4,5; 5,6; 6,7; 7,8; 1,8]; <span class="comment">% Dirichlet bdry edges</span>Neumann   = [];                                       <span class="comment">% Neumann bdry edges</span>mesh = <a href="../../AFEM@matlab/5_Mesh/getmesh.html" class="code" title="function mesh = getmesh(node,elem,Dirichlet,Neumann,type,solu)">getmesh</a>(node,elem,Dirichlet,Neumann);    <span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% PreStep 1: uniform mesh refinement</span><span class="comment">%--------------------------------------------------------------------------</span><span class="keyword">for</span> i=1:2    mesh = <a href="../../AFEM@matlab/4_Refine/uniformrefine.html" class="code" title="function [mesh] = uniformrefine(mesh)">uniformrefine</a>(mesh);<span class="keyword">end</span><span class="comment">% Since our error estimate is ZZ recovery, to be more accurate,</span><span class="comment">% we do not start with very coarse mesh.</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Adaptive Finite Element Method</span><span class="comment">%--------------------------------------------------------------------------</span><span class="keyword">for</span> iter = 1:MaxIt    <span class="comment">% Step1: Solve</span>    [mesh.solu, energy(iter)] = <a href="../../AFEM@matlab/2_Solve/Poisson.html" class="code" title="function [u, energy] = Poisson(mesh, f, g_D, g_N)">Poisson</a>(mesh,@<a href="#_sub1" class="code" title="subfunction z = f(p)">f</a>,@<a href="#_sub2" class="code" title="subfunction z = g_D(p)">g_D</a>,@<a href="#_sub3" class="code" title="subfunction z = g_N(p)">g_N</a>);    view(-47,10);    <span class="comment">% compute the approximation error</span>    l2error(iter) = <a href="../../AFEM@matlab/2_Solve/femerror.html" class="code" title="function [l2error, h1error] = femerror(mesh, u, u_x, u_y)">femerror</a>(mesh,@<a href="#_sub4" class="code" title="subfunction z = u(p)">u</a>,@u_x,@u_y);    pause(0.2)    <span class="comment">% Step2: Estimate</span>    eta = <a href="../../AFEM@matlab/3_Estimate/estimate.html" class="code" title="function eta = estimate(mesh)">estimate</a>(mesh).^2;    <span class="comment">% record number of elements for comparison</span>    N(iter) = length(eta);        <span class="comment">% Step3: Refine</span>    <span class="keyword">for</span> i = 1:RT        [mesh,eta] = <a href="../../AFEM@matlab/4_Refine/bisection.html" class="code" title="function [mesh,eta] = bisection(mesh,eta,theta)">bisection</a>(mesh,eta,theta);    <span class="keyword">end</span>    pause(0.2)        <span class="comment">% Step4: Coarsen</span>    <span class="keyword">if</span> (mod(iter,CF)==0) <span class="comment">% after CF steps of refinement,</span>        <span class="comment">% we do a coarsening to further improve the optimality</span>        [mesh,eta] = <a href="../../AFEM@matlab/4_Refine/coarsening.html" class="code" title="function [mesh,eta] = coarsening(mesh,eta,theta)">coarsening</a>(mesh,eta,theta_c);    <span class="keyword">end</span>        pause(0.2)    <span class="keyword">end</span>[mesh.solu, lastenergy] = <a href="../../AFEM@matlab/2_Solve/Poisson.html" class="code" title="function [u, energy] = Poisson(mesh, f, g_D, g_N)">Poisson</a>(mesh,@<a href="#_sub1" class="code" title="subfunction z = f(p)">f</a>,@<a href="#_sub2" class="code" title="subfunction z = g_D(p)">g_D</a>,@<a href="#_sub3" class="code" title="subfunction z = g_N(p)">g_N</a>);view(-47,10);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Plot convergence rates</span><span class="comment">%--------------------------------------------------------------------------</span>figure(2); set(gcf,<span class="string">'Units'</span>,<span class="string">'normal'</span>); set(gcf,<span class="string">'Position'</span>,[0.02,0.1,0.8,0.6]);subplot(1,2,1)loglog(N,l2error,<span class="string">'-*'</span>, N,N.^(-1),<span class="string">'r-.'</span>); axis tight;title(<span class="string">'L^2 error'</span>, <span class="string">'FontSize'</span>, 14);subplot(1,2,2)loglog(N(1:length(N)-2),energy(1:length(N)-2)-lastenergy,<span class="string">'-*'</span>, <span class="keyword">...</span>       N(1:length(N)-2),N(1:length(N)-2).^(-1),<span class="string">'r-.'</span>); axis tight;title(<span class="string">'Energy error'</span>, <span class="string">'FontSize'</span>, 14);mesh<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% End of function LSHAPE</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Sub functions called by LSHAPE</span><span class="comment">%--------------------------------------------------------------------------</span><a name="_sub1" href="#_subfunctions" class="code">function z = f(p)</a><span class="comment">% load data (right hand side function)</span>z = zeros(size(p,1),1);<span class="comment">%--------------------------------------------------------------------------</span><a name="_sub2" href="#_subfunctions" class="code">function z = g_D(p) </a><span class="comment">% Dirichlet boundary condition</span>z = <a href="#_sub4" class="code" title="subfunction z = u(p)">u</a>(p);<span class="comment">%--------------------------------------------------------------------------</span><a name="_sub3" href="#_subfunctions" class="code">function z = g_N(p) </a><span class="comment">% Neumann boundary condition</span>z = zeros(size(p,1),1);<span class="comment">%--------------------------------------------------------------------------</span><a name="_sub4" href="#_subfunctions" class="code">function z = u(p) </a><span class="comment">% exact solution of the test problem</span>r = sqrt(sum(p.^2,2));theta = atan2(p(:,2),p(:,1));theta = (theta&gt;=0).*theta + (theta&lt;0).*(theta+2*pi);z = r.^(2/3).*sin(2*theta/3);<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> &copy; 2003</address></body></html>
⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -