📄 crack.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 crack</title> <meta name="keywords" content="crack"> <meta name="description" content="CRACK solves Poisson equation in a crack domain with AFEM."> <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 1_Example --><h1>crack</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="box"><strong>CRACK solves Poisson equation in a crack domain with AFEM.</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 crack </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"> CRACK solves Poisson equation in a crack domain with AFEM.</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/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><li><a href="#_sub5" class="code">function z = u_x(p)</a></li><li><a href="#_sub6" class="code">function z = u_y(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 crack</a><span class="comment">% CRACK solves Poisson equation in a crack domain with AFEM.</span><span class="comment">%</span><span class="comment">% L. Chen & 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.4; theta_c = 0.1; <span class="comment">% theta and theta_c are parameters for bisection and coarsening</span>RT = 1; CF = 4; 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 = sparse(MaxIt,1); energy = sparse(MaxIt,1); l2error = sparse(MaxIt,1); h1error = sparse(MaxIt,1);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% PreStep 0: generate initial mesh</span><span class="comment">%--------------------------------------------------------------------------</span>node = [1,0; 0,1; -1,0; 0,-1; 0,0; 1,0]; <span class="comment">% nodes</span>elem = [5,1,2; 5,2,3; 5,3,4; 5,4,6]; <span class="comment">% elements</span>Dirichlet = [1,2; 2,3; 3,4; 4,6; 1,5; 5,6]; <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 refinement</span><span class="comment">%--------------------------------------------------------------------------</span><span class="keyword">for</span> i=1:3 mesh = <a href="../../AFEM@matlab/4_Refine/bisection.html" class="code" title="function [mesh,eta] = bisection(mesh,eta,theta)">bisection</a>(mesh,ones(size(mesh.elem,1),1),1);<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>); <span class="comment">% compute error for comparison</span> [l2error(iter),h1error(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>,@<a href="#_sub5" class="code" title="subfunction z = u_x(p)">u_x</a>,@<a href="#_sub6" class="code" title="subfunction z = u_y(p)">u_y</a>); 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 <span class="comment">% the inner loop will refine the mesh without solving</span> [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><span class="comment">% Solve the equation on the finest mesh</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>);<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,3,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,3,2)loglog(N,h1error,<span class="string">'-*'</span>, N,N.^(-.5),<span class="string">'r-.'</span>); axis tight;title(<span class="string">'H^1 error'</span>, <span class="string">'FontSize'</span>, 14);subplot(1,3,3)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 CRACK</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Sub functions called by CRACK</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 = ones(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</span>r = sqrt(sum(p.^2,2));z = sqrt(0.5*(r-p(:,1)))-0.25*r.^2;<span class="comment">%--------------------------------------------------------------------------</span><a name="_sub5" href="#_subfunctions" class="code">function z = u_x(p) </a><span class="comment">% x-derivative of the exact solution</span>r = sqrt(sum(p.^2,2));z = (p(:,1)./r-1)./sqrt(8*(r-p(:,1)))-0.5*p(:,1);<span class="comment">%--------------------------------------------------------------------------</span><a name="_sub6" href="#_subfunctions" class="code">function z = u_y(p) </a><span class="comment">% y-derivative of the exact solution</span>r = sqrt(sum(p.^2,2));z = p(:,2)./r./sqrt(8*(r-p(:,1)))-0.5*p(:,2);<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 + -