fem_1st_to_higher_order

 FEM_1ST_TO_HIGH_ORDER:  Modify the FEM for high order FEM called as
    [bound,elem,nodes] = fem_1st_to_higher_order( fwd_model )
 where field fwd_model.approx_type = 'tri3'  - 2D linear
                                   = 'tri6'  - 2D quadratic
                                   = 'tri10' - 2D cubic
                                   = 'tet4'  - 3D linear
                                   = 'tet10' - 3D quadratic
 fwd_model : is a fwd_model structure
 boundary/elems : new boundary and element numbers in connectivity matrix
 nodes : fem nodes including extra nodes added in refinement

M Crabb - 29.06.2012

0001 function [boundary,elems,nodes]=fem_1st_to_higher_order(fwd_model)
0002 % FEM_1ST_TO_HIGH_ORDER:  Modify the FEM for high order FEM called as
0003 %    [bound,elem,nodes] = fem_1st_to_higher_order( fwd_model )
0004 % where field fwd_model.approx_type = 'tri3'  - 2D linear
0005 %                                   = 'tri6'  - 2D quadratic
0006 %                                   = 'tri10' - 2D cubic
0007 %                                   = 'tet4'  - 3D linear
0008 %                                   = 'tet10' - 3D quadratic
0009 % fwd_model : is a fwd_model structure
0010 % boundary/elems : new boundary and element numbers in connectivity matrix
0011 % nodes : fem nodes including extra nodes added in refinement
0012 %
0013 %M Crabb - 29.06.2012
0014 
0015 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return ; end
0016 
0017 copt.cache_obj = {fwd_model.nodes, fwd_model.elems, ...
0018              fwd_model.approx_type, fwd_model.boundary};
0019 copt.fstr = 'boundaryN_elemsN_nodes';
0020 
0021 [boundary,elems,nodes]= eidors_cache(@mc_fem_modify,fwd_model,copt);
0022 end
0023 
0024 function [boundary,elems,nodes]=mc_fem_modify(mdl)
0025 %Function  that takes a fwd_model structure (with at least fwd_model.boundary,
0026 %fwd_model.elems, fwd_model.nodes), reorders these to form anti-clockwise sets, and
0027 %then adds extra nodes for p-refinement.
0028 
0029 %Reorder element nodes anti-clockwise relative to first node
0030 mdl = linear_elem_reorder(mdl,-1);
0031 
0032 %Reorder boundary nodes anti-clockwise relative to interior node
0033 mdl = linear_bound_reorder(mdl,-1);
0034 
0035 %Find nodes, elems, boundary and problem dimensions
0036 n_nodes = size(mdl.nodes,1); n_elems = size(mdl.elems,1); n_bounds=size(mdl.boundary,1);
0037 nodedim = size(mdl.nodes,2);
0038 
0039 %Copy nodes, elems and boundaries
0040 nodes_n=mdl.nodes; boundary_n=mdl.boundary; elems_n=mdl.elems;
0041 
0042 %Compare strings
0043 if(strcmp(mdl.approx_type,'tri6'))    
0044   %No nodes (really vertices) per element and boundary
0045   n_nodes_per_elem = 3; 
0046   n_nodes_per_boundary=1;
0047 elseif(strcmp(mdl.approx_type,'tri10'))
0048   %No nodes (really vertices) per element and boundary
0049   n_nodes_per_elem = 7; 
0050   n_nodes_per_boundary=2;
0051 elseif(strcmp(mdl.approx_type,'tet10'));
0052   %No nodes (really vertices) per element and boundary
0053   n_nodes_per_elem = 6; 
0054   n_nodes_per_boundary=3;
0055 else
0056     error('mc_fem_modify: Element type ("%s") not recognised',mdl.approx_type);
0057 end
0058 
0059 %Maximum number of new nodes i.e. unconnected elemets and boundaries
0060 n_nodes_e = n_nodes + n_nodes_per_elem*n_elems;
0061 n_nodes_e_b = n_nodes_e + n_nodes_per_boundary*n_bounds;
0062 
0063 %Now reshape the nodes, elems and boundary for extra FEM nodes
0064 nodes_n(n_nodes+1:n_nodes_e_b,:) = 0;
0065 elems_n(:,nodedim+2:nodedim+1+n_nodes_per_elem) = reshape((n_nodes+1):n_nodes_e,n_nodes_per_elem,n_elems)';
0066 boundary_n(:,nodedim+1:nodedim+n_nodes_per_boundary) = reshape((n_nodes_e+1):n_nodes_e_b,n_nodes_per_boundary,n_bounds)';
0067 
0068 %Calculate positions of now nodes dependent on approx_type
0069 if(strcmp(mdl.approx_type,'tri6'))    
0070   nodes_n(elems_n(:,4),:)  = 0.5*(nodes_n(elems_n(:,1),:) + nodes_n(elems_n(:,2),:));
0071   nodes_n(elems_n(:,5),:)  = 0.5*(nodes_n(elems_n(:,2),:) + nodes_n(elems_n(:,3),:));
0072   nodes_n(elems_n(:,6),:)  = 0.5*(nodes_n(elems_n(:,3),:) + nodes_n(elems_n(:,1),:));
0073   nodes_n(boundary_n(:,3),:) = 0.5*(nodes_n(boundary_n(:,1),:) + nodes_n(boundary_n(:,2),:));
0074 elseif(strcmp(mdl.approx_type,'tri10'))
0075   nodes_n(elems_n(:,4),:) = 2.0/3.0*nodes_n(elems_n(:,1),:) + 1.0/3.0*nodes_n(elems_n(:,2),:);
0076   nodes_n(elems_n(:,5),:) = 1.0/3.0*nodes_n(elems_n(:,1),:) + 2.0/3.0*nodes_n(elems_n(:,2),:);
0077   nodes_n(elems_n(:,6),:) = 2.0/3.0*nodes_n(elems_n(:,2),:) + 1.0/3.0*nodes_n(elems_n(:,3),:);
0078   nodes_n(elems_n(:,7),:) = 1.0/3.0*nodes_n(elems_n(:,2),:) + 2.0/3.0*nodes_n(elems_n(:,3),:);
0079   nodes_n(elems_n(:,8),:) = 2.0/3.0*nodes_n(elems_n(:,3),:) + 1.0/3.0*nodes_n(elems_n(:,1),:);
0080   nodes_n(elems_n(:,9),:) = 1.0/3.0*nodes_n(elems_n(:,3),:) + 2.0/3.0*nodes_n(elems_n(:,1),:);
0081   nodes_n(elems_n(:,10),:) = 1.0/3.0*nodes_n(elems_n(:,1),:) + 1.0/3.0*nodes_n(elems_n(:,2),:) + 1.0/3.0*nodes_n(elems_n(:,3),:);
0082   nodes_n(boundary_n(:,3),:) = 2.0/3.0*nodes_n(boundary_n(:,1),:) + 1.0/3.0*nodes_n(boundary_n(:,2),:);
0083   nodes_n(boundary_n(:,4),:) = 1.0/3.0*nodes_n(boundary_n(:,1),:) + 2.0/3.0*nodes_n(boundary_n(:,2),:);
0084 elseif(strcmp(mdl.approx_type,'tet10'))
0085   nodes_n(elems_n(:,5),:) = 0.5*nodes_n(elems_n(:,1),:) + 0.5*nodes_n(elems_n(:,2),:);
0086   nodes_n(elems_n(:,6),:) = 0.5*nodes_n(elems_n(:,1),:) + 0.5*nodes_n(elems_n(:,3),:);
0087   nodes_n(elems_n(:,7),:) = 0.5*nodes_n(elems_n(:,1),:) + 0.5*nodes_n(elems_n(:,4),:);
0088   nodes_n(elems_n(:,8),:) = 0.5*nodes_n(elems_n(:,2),:) + 0.5*nodes_n(elems_n(:,3),:);
0089   nodes_n(elems_n(:,9),:) = 0.5*nodes_n(elems_n(:,3),:) + 0.5*nodes_n(elems_n(:,4),:);
0090   nodes_n(elems_n(:,10),:) = 0.5*nodes_n(elems_n(:,2),:) + 0.5*nodes_n(elems_n(:,4),:);
0091   nodes_n(boundary_n(:,4),:) = 0.5*nodes_n(boundary_n(:,1),:) + 0.5*nodes_n(boundary_n(:,2),:);
0092   nodes_n(boundary_n(:,5),:) = 0.5*nodes_n(boundary_n(:,1),:) + 0.5*nodes_n(boundary_n(:,3),:);
0093   nodes_n(boundary_n(:,6),:) = 0.5*nodes_n(boundary_n(:,2),:) + 0.5*nodes_n(boundary_n(:,3),:);
0094 else
0095     error('mc_fem_modify: Element type ("%s") not recognised',mdl.approx_type);
0096 end
0097 
0098 %Find the unique nodes by row
0099 [nodes, a, b] = unique(nodes_n(n_nodes+1:end,:),'rows');
0100 
0101 %Add the unique nodes to the model
0102 nodes_n = [mdl.nodes; nodes];
0103 
0104 %Now find unique nodes
0105 c = [1:n_nodes n_nodes+b']; 
0106 elems_n = c(elems_n);
0107 boundary_n = c(boundary_n);
0108 
0109 %Reassign matrices
0110 nodes=nodes_n; elems=elems_n; boundary=boundary_n;
0111 
0112 end
0113 
0114 %Function to reorder nodes consistently within element/boundary
0115 %Elem Default : Arrange nodes locally in mdl.elems so that:
0116 %2D: Counter-clockwise looking down onto element
0117 %3D: Counter-clockwise looking into element from first node
0118 %    (clockwise looking down from ouside element opposite first node)
0119 %Boundary Default: Arrange nodes locally in mdl.boundary so that:
0120 %2D: Counter-clockwise looking down relative to interior node
0121 %3D: Clockwise relative looking into element from interior node
0122 %    (counter-clockwise looking down from outside element opposite interior node)
0123 %Conventions MUST match reference element geometry in system_mat_higher_order
0124 function [mdl] = linear_elem_reorder(mdl,ccw)
0125     %Find elements, nodes and no. of elements
0126     elemstruc=mdl.elems;  nodestruc=mdl.nodes; nelems=size(elemstruc,1);
0127     for e=1:nelems;
0128         %Row vector of the node numbers and the no. of vertices
0129         enodes = elemstruc(e,:); elenode = size(enodes,2);
0130     
0131         %Matrix of nodal positions [elenodexdim] (Linear dimension==elenode-1)
0132         nd = nodestruc(enodes,:);
0133     
0134        %Calculate area(2D)/volume(3D) defined by the vertices
0135         area= det([ones(elenode,1),nd]); areasign=sign(area);
0136     
0137         %If sign is (pos) neg swap two nodes (last two will suffice..)
0138         if(areasign == ccw) %Swap last two entries of enodes
0139             temp=enodes(elenode-1);
0140             enodes(elenode-1)=enodes(elenode);
0141             enodes(elenode) = temp;
0142         end
0143         elemstruc(e,:)=enodes; %Put enodes back into elementnodes matrix
0144     end
0145     %Reassign the elements
0146     mdl.elems=elemstruc;
0147 end
0148 
0149 function [mdl]=linear_bound_reorder(mdl,ccw)
0150     %Reorder nodes on each electrode's boundary
0151     %SPEED UP
0152     %1fix_model with option opt.bound_elec2elem=1
0153     
0154     %Find boundary, elems, nodes, elecs
0155     boundstruc=mdl.boundary; elemstruc=mdl.elems; nodestruc=mdl.nodes; elecstruc=mdl.electrode;
0156     %Find no. elecs and node dimension
0157     nelecs=size(elecstruc,2); nodedim=size(nodestruc,2);
0158 
0159     %Reorder boundaries belonging to electrodes
0160     for ke=1:nelecs
0161         %The boundary numbers and areas, outputs rows of mdl.boundary of electrode
0162         bdy_idx=find_electrode_bdy(boundstruc(:,1:nodedim),nodestruc,elecstruc(ke).nodes);                 
0163             
0164         for ii=1:length(bdy_idx);  
0165             %Get bdy idx
0166             bb=bdy_idx(ii);
0167                                     
0168             %Vector of vertes numbers of this boundary
0169             bbnodes=boundstruc(bb,:);
0170             if(nodedim==2) %2D problem
0171                 %Find row(s) of elems for which each boundary vertex belongs
0172                 [rownode1,blah]=find(elemstruc==bbnodes(1));
0173                 [rownode2,blah]=find(elemstruc==bbnodes(2));
0174         
0175                 %Intersection of rownode1 and rownode2 is the unique element
0176                 boundiielem=intersect(rownode1,rownode2);
0177             elseif(nodedim==3) %3D problem
0178                 %Find row(s) of elems for which each node belongs
0179                 [rownode1,blah]=find(elemstruc==bbnodes(1));
0180                 [rownode2,blah]=find(elemstruc==bbnodes(2));
0181                 [rownode3,blah]=find(elemstruc==bbnodes(3));
0182 
0183                 %Intersection of rownode1 and rownode2 gives a choic
0184                 rownode1node2=intersect(rownode1,rownode2);
0185 
0186                 %Intersection of rownode3 with vector above is unique element
0187                 boundiielem=intersect(rownode3,rownode1node2);  
0188             end
0189             %Store this unique number in a vecto
0190             %FIXME!!! 3D Inclusion models seem to include the boundaries of the
0191             %inclusions as boundaries. In this case there may be multiple
0192             %boundaries. Just choose first one for time being, and if this
0193             %is ACTUAL boundary, then should be unique.....
0194             beleind=boundiielem(1);     
0195             
0196             %Coordinate of nodes of boundaries element
0197             belenodes=elemstruc(beleind,:);
0198   
0199             %Find unique internal node and order element node so internal is first
0200             intnode=setdiff(belenodes,bbnodes); elemboundnode=[intnode,bbnodes];
0201     
0202             %List by row coordinates of the element
0203             nodepos=nodestruc(elemboundnode,:); nvertices=size(belenodes,2);     
0204 
0205             %Calculate area(2D)/volume(3D) defined by the vertices
0206             area= det([ones(nvertices,1),nodepos]); areasign=sign(area);
0207     
0208             %If positive area, swap the last two nodes
0209             if(areasign == -ccw) %Swap last two entries of enodes
0210                 temp=elemboundnode(end-1);
0211                 elemboundnode(end-1)=elemboundnode(end);
0212                 elemboundnode(end) = temp;
0213             end
0214     
0215             %Put the node numbers back into mdl.bound(bb) (excluding internal node)
0216             boundstruc(bb,:)=elemboundnode(2:end);
0217         end              
0218     end
0219     %Reassign the boundary
0220     mdl.boundary=boundstruc;
0221 end
0222 
0223 function do_unit_test
0224      do_unit_test_2D;
0225      do_unit_test_3D;
0226 end
0227 
0228 function do_unit_test_2D
0229     imdl=mk_common_model('c2C',16);
0230     fmdl=imdl.fwd_model;
0231     fmdl.approx_type='tri6';
0232     [bou,ele,nod]=fem_1st_to_higher_order(fmdl);
0233 
0234     unit_test_cmp('NOD',nod(1:5,:), ...
0235 [                  0                   0
0236   -0.005450260769179   0.083154910269884
0237    0.083154910269884   0.005450260769179
0238    0.005450260769179  -0.083154910269884
0239   -0.083154910269884  -0.005450260769179],1e-8);
0240     unit_test_cmp('ELE',ele(1:5,1:6), ...
0241 [1, 3, 2, 783, 780, 755
0242  1, 4, 3, 760, 785, 783
0243  1, 5, 4, 732, 735, 760
0244  1, 2, 5, 755, 730, 732
0245  7, 2, 3, 792, 780, 810]);
0246 
0247 end
0248 
0249 function do_unit_test_3D
0250     fmdl=mk_library_model('adult_male_16el');
0251     fmdl.stimulation = mk_stim_patterns(16,1,'{ad}','{ad}');
0252     fmdl.approx_type='tet10';
0253     [bou,ele,nod]=fem_1st_to_higher_order(fmdl);
0254 
0255     unit_test_cmp('NOD',nod(1:5,:), ...
0256 [ -0.812999999999999, -0.439800000000001, 0
0257   -0.720700000000001, -0.494600000000000, 0
0258   -0.886900000000000, -0.361400000000000, 0
0259   -0.616600000000001, -0.523199999999999, 0
0260   -0.516400000000002, -0.562899999999999, 0],1e-8 );
0261 
0262 end