h_refine_square_domain

Perform uniform h refinement of square CEM

0001 function mdl = h_refine(mdl)
0002 %Perform uniform h refinement of square CEM
0003 
0004 %Copy nodes, elems and boundaries
0005 nodes=mdl.nodes; boundary=mdl.boundary; elems=mdl.elems;
0006 
0007 %Find nodes, elems, boundary and problem dimensions
0008 n_nodes = size(nodes,1); n_elems = size(elems,1); n_boundary=size(boundary,1);
0009 
0010 %Maximum number of new nodes i.e. unconnected elemets and boundaries
0011 n_nodes_new = n_nodes + 3*n_elems + n_boundary; %Maximum number of new nodes
0012 n_boundary_new = n_boundary*2; %Number of new boundaries EXACT
0013 n_elems_new = n_elems*4; %Number of new elements EXACT
0014 
0015 %Storage for these new arrays
0016 nodes_new=zeros(n_nodes_new,2);
0017 elems_new=zeros(n_elems_new,3);
0018 boundary_new=zeros(n_boundary_new,2);
0019 
0020 %Keep original nodes in n_nodes_new
0021 nodes_new(1:n_nodes,:)=nodes;
0022 
0023 %Loop through each element and fine edge midpoints.
0024 cnt=0;
0025 for i=1:n_elems
0026 %Find new nodes on midpoints
0027 nodes_new(n_nodes+(i-1)*3+1,:) = 0.5*(nodes(elems(i,1),:) + nodes(elems(i,2),:));
0028 nodes_new(n_nodes+(i-1)*3+2,:) = 0.5*(nodes(elems(i,1),:) + nodes(elems(i,3),:));
0029 nodes_new(n_nodes+(i-1)*3+3,:) = 0.5*(nodes(elems(i,2),:) + nodes(elems(i,3),:));
0030 
0031 %Create the four new elements
0032 elems_new((i-1)*4+1,1) = elems(i,1); elems_new((i-1)*4+1,2) = n_nodes+3*(i-1)+1; elems_new((i-1)*4+1,3)=n_nodes+3*(i-1)+2;
0033 elems_new((i-1)*4+2,1) = elems(i,2); elems_new((i-1)*4+2,2) = n_nodes+3*(i-1)+1; elems_new((i-1)*4+2,3)=n_nodes+3*(i-1)+3;
0034 elems_new((i-1)*4+3,1) = elems(i,3); elems_new((i-1)*4+3,2) = n_nodes+3*(i-1)+2; elems_new((i-1)*4+3,3)=n_nodes+3*(i-1)+3;
0035 elems_new((i-1)*4+4,1) = n_nodes+3*(i-1)+1; elems_new((i-1)*4+4,2) = n_nodes+3*(i-1)+2; elems_new((i-1)*4+4,3)= n_nodes+3*(i-1)+3;
0036 end
0037 
0038 %Loop through each boundary and fine edge midpoints
0039 for i=1:n_boundary
0040 %Find node and add a row of
0041 nodes_new(n_nodes+3*n_elems+i,:) = 0.5*(nodes(boundary(i,1),:) + nodes(boundary(i,2),:));
0042 
0043 %Create the two new boundarys
0044 boundary_new((i-1)*2+1,1) = boundary(i,1); boundary_new((i-1)*2+1,2) = n_nodes+3*n_elems+i; 
0045 boundary_new((i-1)*2+2,1) = n_nodes+3*n_elems+i; boundary_new((i-1)*2+2,2) = boundary(i,2); 
0046 end
0047 
0048 %Find the unique nodes by row
0049 [nodes, a, b] = unique(nodes_new(n_nodes+1:end,:),'rows');
0050 
0051 %Add the unique nodes to the model
0052 nodes_new = [mdl.nodes; nodes];
0053 
0054 %Now find unique nodes
0055 c = [1:n_nodes n_nodes+b']; 
0056 elems_new = c(elems_new);
0057 boundary_new = c(boundary_new);
0058 
0059 %Reassign matrices
0060 nodes=nodes_new; elems=elems_new; boundary=boundary_new;
0061 
0062 %Find new CEM electrode nodes
0063 %Recalculate the number of nodes
0064 n_nodes = size(nodes,1);
0065 cntE1 = 0; cntE2 = 0;
0066 %Loop through boundary and see if electrode is in between electrodes
0067 for i=1:n_nodes
0068     if(nodes(i,2)==0) %We are on bottom boundary
0069         if( pi/5 <= nodes(i,1) && nodes(i,1)  <= 2*pi/5)
0070             cntE1 = cntE1+1;
0071             E1node(cntE1) = i;
0072         elseif(  3*pi/5 <= nodes(i,1) && nodes(i,1)  <= 4*pi/5)
0073             cntE2 = cntE2+1;
0074             E2node(cntE2) = i;                
0075         end
0076     end            
0077 end
0078 %Electrode
0079 mdl.electrode(1).nodes=E1node;
0080 mdl.electrode(2).nodes=E2node;        
0081 
0082 %Reassign
0083 mdl.nodes=nodes; mdl.elems=elems; mdl.boundary=boundary;
0084     
0085 end