📄 bisection.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 bisection</title> <meta name="keywords" content="bisection"> <meta name="description" content="BISECTION refines the triangulation using newest vertex bisection"> <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 4_Refine --><h1>bisection</h1><h2><a name="_name"></a>PURPOSE <a href="#_top"><img alt="^" border="0" src="../../up.png"></a></h2><div class="box"><strong>BISECTION refines the triangulation using newest vertex bisection</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 [mesh,eta] = bisection(mesh,eta,theta) </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"> BISECTION refines the triangulation using newest vertex bisection USAGE [mesh,eta] = bisection(mesh,eta,theta) INPUT mesh: current mesh eta: error indicator for each triangle theta: parameter in (0,1). We mark minimal number of triangles M such that \sum_{T \in M} \eta_T > \theta*\sum\eta_T OUTPUT mesh: new mesh after refinement eta: new error indicator for each triangle REFERENCE Long Chen, Short bisection implementation in MATLAB Research Notes, 2006</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/Lbig.html" class="code" title="function Lbig">Lbig</a> LBIG solves Poisson equation in a L-shaped domain with FEM with</li><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/circle.html" class="code" title="function circle">circle</a> CIRCLE simulates adaptive grids for a problem with moving circlular</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></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 bdEdge = updatebd(bdEdge,marker,d2p)</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 [mesh,eta] = bisection(mesh,eta,theta)</a><span class="comment">% BISECTION refines the triangulation using newest vertex bisection</span><span class="comment">%</span><span class="comment">% USAGE</span><span class="comment">% [mesh,eta] = bisection(mesh,eta,theta)</span><span class="comment">%</span><span class="comment">% INPUT</span><span class="comment">% mesh: current mesh</span><span class="comment">% eta: error indicator for each triangle</span><span class="comment">% theta: parameter in (0,1).</span><span class="comment">% We mark minimal number of triangles M such that</span><span class="comment">% \sum_{T \in M} \eta_T > \theta*\sum\eta_T</span><span class="comment">%</span><span class="comment">% OUTPUT</span><span class="comment">% mesh: new mesh after refinement</span><span class="comment">% eta: new error indicator for each triangle</span><span class="comment">%</span><span class="comment">% REFERENCE</span><span class="comment">% Long Chen,</span><span class="comment">% Short bisection implementation in MATLAB</span><span class="comment">% Research Notes, 2006</span><span class="comment">%</span><span class="comment">% L. Chen & C. Zhang 11-12-2006</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Construct data structure</span><span class="comment">%--------------------------------------------------------------------------</span>edge = [mesh.elem(:,[1,2]); mesh.elem(:,[1,3]); mesh.elem(:,[2,3])];edge = unique(sort(edge,2),<span class="string">'rows'</span>);N = size(mesh.node,1);NT = size(mesh.elem,1); NE = size(edge,1);dualEdge = sparse(mesh.elem(:,[1,2,3]),mesh.elem(:,[2,3,1]),[1:NT,1:NT,1:NT]);d2p = sparse(edge(:,[1,2]),edge(:,[2,1]),[1:NE,1:NE]);<span class="comment">% Detailed explanation can be founded at</span><span class="comment">% Manual --> Data Structure --> Auxlliary data structure</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Meomery management for node arrary</span><span class="comment">%--------------------------------------------------------------------------</span>recycle = find(mesh.type==0); last = length(recycle);<span class="comment">% Coarsening can create empty spaces in the node array. We collect those</span><span class="comment">% scattered spaces in recycle arrary and 'last' will point out the last</span><span class="comment">% empty node index.</span><span class="comment">% mesh.type array is used to distinguish the type of nodes:</span><span class="comment">% 0: empty nodes (deleted by coarsening);</span><span class="comment">% 1: nodes in the initial triangulation or nodes from regular refinement;</span><span class="comment">% 2: new added nodes due to refinement;</span><span class="comment">% 5: temporay node (will be deleted when bisection finished).</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Mark triangles according to the error indicator</span><span class="comment">%--------------------------------------------------------------------------</span>total = sum(eta); current = 0; [temp,ix] = sort(-eta); <span class="comment">% sort in descent order</span>marker = zeros(NE,1); <span class="comment">% initialize for possible new nodes</span>mesh.node = [mesh.node; zeros(NE,2)];mesh.type = [mesh.type; uint8(5*ones(NE,1))];mesh.solu = [mesh.solu; zeros(NE,1)];<span class="keyword">for</span> t = 1:NT <span class="keyword">if</span> (current > theta*total), <span class="keyword">break</span>; <span class="keyword">end</span> <span class="comment">% est on marked elem big enough</span> index = 1; ct = ix(t); <span class="keyword">while</span> (index==1) base = d2p(mesh.elem(ct,2),mesh.elem(ct,3)); <span class="keyword">if</span> (marker(base)>0) <span class="comment">% base is already marked</span> index = 0; <span class="keyword">else</span> current = current + eta(ct); <span class="keyword">if</span> (last==0) newNode = N + 1; N = N+1; <span class="keyword">end</span> <span class="keyword">if</span> (last>0) newNode = recycle(last); last = last-1; <span class="keyword">end</span> marker( d2p(mesh.elem(ct,2),mesh.elem(ct,3)) ) = newNode; <span class="comment">% A new node is added to the mesh. Numerical solution at this</span> <span class="comment">% new added node is approximated by linear interpolation. Type</span> <span class="comment">% of this node is 2 (generated by newest vertex bisection)</span> mesh.node(newNode,:) = ( mesh.node(mesh.elem(ct,2),:) + <span class="keyword">...</span> mesh.node(mesh.elem(ct,3),:) )/2; mesh.solu(newNode) = ( mesh.solu(mesh.elem(ct,2)) + <span class="keyword">...</span> mesh.solu(mesh.elem(ct,3)) )/2; mesh.type(newNode) = 2; <span class="comment">% Find the element which shares the base edge of the current</span> <span class="comment">% element. If it is 0, it means the base of the current element</span> <span class="comment">% is on the boundary.</span> ct = dualEdge( mesh.elem(ct,3), mesh.elem(ct,2) ); <span class="keyword">if</span> ( ct == 0 ), index = 0; <span class="keyword">end</span> <span class="comment">% the while will ended if</span> <span class="comment">% 1. ct==0 means we are on the boundary</span> <span class="comment">% 2. base(ct) is already marked</span> <span class="keyword">end</span> <span class="keyword">end</span> <span class="comment">% end while for recursive marking</span><span class="keyword">end</span> <span class="comment">% end of for loop on all elements</span><span class="comment">% Detailed explanation of the algorithm can be found at</span><span class="comment">% Manual --> Algorithms --> Bisection</span><span class="comment">% delete possible empty entries</span>ix = (mesh.type == 5); mesh.node(ix,:) = []; mesh.type(ix) = [];mesh.solu(ix) = [];<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Refine marked edges for each triangle</span><span class="comment">%--------------------------------------------------------------------------</span>numnew = 2*sum(marker~=0); <span class="comment">% number of new elements need to be added</span>mesh.elem = [mesh.elem; zeros(numnew,3)];eta = [eta; zeros(numnew,1)];inew = NT + 1; <span class="comment">% index for current new added right child</span><span class="keyword">for</span> t = 1:NT base = d2p(mesh.elem(t,2),mesh.elem(t,3)); <span class="keyword">if</span> (marker(base)>0) p = [mesh.elem(t,:), marker(base)]; <span class="comment">% Case 1: divide the current marked triangle</span> mesh.elem(t,:) = [p(4),p(1),p(2)]; <span class="comment">% t is always a left child</span> mesh.elem(inew,:) = [p(4),p(3),p(1)]; <span class="comment">% new is a right child</span> eta(t) = eta(t)/2; <span class="comment">% update error indicators</span> eta(inew) = eta(t); inew = inew + 1; <span class="comment">% Case 2: divide the right child, different, careful!!!</span> right = d2p(p(3),p(1)); <span class="keyword">if</span> (marker(right)>0) mesh.elem(inew-1,:) = [marker(right),p(4),p(3)]; mesh.elem(inew,:) = [marker(right),p(1),p(4)]; eta(inew-1) = eta(inew-1)/2; eta(inew) = eta(inew-1); inew = inew + 1; <span class="keyword">end</span> <span class="comment">% Case 3: divide the left child, similar to the case 1.</span> left = d2p(p(1),p(2)); <span class="keyword">if</span> (marker(left)>0) mesh.elem(t,:) = [marker(left),p(4),p(1)]; mesh.elem(inew,:) = [marker(left),p(2),p(4)]; eta(t) = eta(t)/2; eta(inew) = eta(t); inew = inew + 1; <span class="keyword">end</span> <span class="keyword">end</span> <span class="comment">% end of refinement of one element</span><span class="keyword">end</span> <span class="comment">% end of for loop on all elements</span><span class="comment">% delete possible empty entries</span>mesh.elem = mesh.elem(1:inew-1,:);eta = eta(1:inew-1,:);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Update boundary edges</span><span class="comment">%--------------------------------------------------------------------------</span>mesh.Dirichlet = <a href="#_sub1" class="code" title="subfunction bdEdge = updatebd(bdEdge,marker,d2p)">updatebd</a>(mesh.Dirichlet,marker,d2p);mesh.Neumann = <a href="#_sub1" class="code" title="subfunction bdEdge = updatebd(bdEdge,marker,d2p)">updatebd</a>(mesh.Neumann,marker,d2p);<span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% End of function BISECTION</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% Sub functions called by BISECTION</span><span class="comment">%--------------------------------------------------------------------------</span><a name="_sub1" href="#_subfunctions" class="code">function bdEdge = updatebd(bdEdge,marker,d2p)</a><span class="comment">% UPDATEDBD refine the boundary edges</span><span class="comment">%</span><span class="comment">% USAGE</span><span class="comment">% bdEdge = updatebd(bdEdge,marker,d2p)</span><span class="comment">%</span><span class="comment">% INPUT</span><span class="comment">% bdEdge: set of boundary edges</span><span class="comment">% marker: new node index for marked edge</span><span class="comment">% d2p: index mapping from dual edge to primary edge</span><span class="comment">%</span><span class="comment">% OUTPUT</span><span class="comment">% bdEdge: set of refined boundary edges</span><span class="comment">%</span>NB = size(bdEdge,1);<span class="keyword">for</span> k = 1:NB i = bdEdge(k,1); j = bdEdge(k,2); <span class="keyword">if</span> marker(d2p(i,j)) >0 bdEdge(k,:) = [i,marker(d2p(i,j))]; bdEdge(size(bdEdge,1)+1,:) = [marker(d2p(i,j)),j]; <span class="keyword">end</span><span class="keyword">end</span><span class="comment">%--------------------------------------------------------------------------</span><span class="comment">% End of function UPDATEDBD</span><span class="comment">%--------------------------------------------------------------------------</span></pre></div><hr><address>Generated on Fri 17-Nov-2006 11:02:53 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 + -