fix_model

 FIX_MODEL: Add useful fields to a model
    [mdl] = fix_model(mdl,options)
 INPUT:
    mdl - an FEM model with at least the following fields:
       .name
       .nodes
       .elem
    options - not coded yet

 OUTPUT:
    mdl - a copy of the input model with these additional fields:
       .boundary
       .boundary_face
       .faces
       .face2elem
       .elem2face
       .elem_centre
       .face_centre
       .normals
       .max_edge_len (per elem)

 The elems will be reordered so that all faces are counter-clockwise.

0001 function [mdl] = fix_model(mdl,options)
0002 % FIX_MODEL: Add useful fields to a model
0003 %    [mdl] = fix_model(mdl,options)
0004 % INPUT:
0005 %    mdl - an FEM model with at least the following fields:
0006 %       .name
0007 %       .nodes
0008 %       .elem
0009 %    options - not coded yet
0010 %
0011 % OUTPUT:
0012 %    mdl - a copy of the input model with these additional fields:
0013 %       .boundary
0014 %       .boundary_face
0015 %       .faces
0016 %       .face2elem
0017 %       .elem2face
0018 %       .elem_centre
0019 %       .face_centre
0020 %       .normals
0021 %       .max_edge_len (per elem)
0022 %
0023 % The elems will be reordered so that all faces are counter-clockwise.
0024 
0025 % (C) 2011 Bartlomiej Grychtol. Licensed under GPL v2 or v3
0026 % $Id: fix_model.html 2819 2011-09-07 16:43:11Z aadler $
0027 
0028 if isstr(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0029 
0030 
0031 mdl = linear_reorder(mdl); %counter-clockwise
0032 if ~isfield(mdl,'boundary')
0033     mdl.boundary = find_boundary(mdl);
0034 end
0035 [mdl.faces mdl.face2elem mdl.elem2face] = calc_faces(mdl);
0036 mdl.boundary_face = mdl.face2elem(:,2)==0;
0037 mdl.elem_centre = interp_mesh(mdl, 0);
0038 tmp = mdl;
0039 tmp.elems = tmp.faces;
0040 mdl.face_centre = interp_mesh(tmp,0);
0041 mdl.normals = calc_normals(mdl);
0042 mdl.inner_normal = test_inner_normal( mdl );
0043 mdl.max_edge_len = calc_longest_edge(mdl.elems,mdl.nodes);
0044 mdl.elem_volume = get_elem_volume(mdl);
0045 
0046 
0047 % decrease memory footprint
0048 mdl.elems = uint32(mdl.elems);
0049 mdl.faces = uint32(mdl.faces);
0050 mdl.elem2face = uint32(mdl.elem2face);
0051 mdl.face2elem = uint32(mdl.face2elem);
0052 
0053 % Test whether normal points into or outsize
0054 % mdl.inner_normal(i,j) = 1 if face i of elem j points in
0055 function inner_normal = test_inner_normal( mdl );
0056    inner_normal = logical(zeros(size(mdl.elem2face)));
0057    d = elem_dim(mdl) + 1;
0058    for i=1:num_elems(mdl);
0059       el_faces = mdl.elem2face(i,:);
0060       el_ctr   = repmat( mdl.elem_centre(i,:), d, 1);
0061       vec_fa_el= el_ctr -  mdl.face_centre(el_faces,:);
0062       normal_i  = mdl.normals(el_faces,:);
0063       dot_prod = sum( normal_i.*vec_fa_el, 2 );
0064       inner_normal(i,:) = dot_prod' > 0;
0065    end
0066 
0067 function [faces, face2elem, elem2face] = calc_faces(mdl)
0068 
0069 
0070 e_dim = mdl_dim(mdl);
0071 
0072 
0073 %elem2face = zeros(size(mdl.elems));
0074 switch e_dim
0075     case 2
0076         idx = [1 2; 2 3; 1 3];
0077     case 3
0078         idx = [1 2 3; 2 3 4; 1 2 4; 1 3 4];
0079 end
0080 elem_sorted = sort(mdl.elems,2);
0081 [faces ib ia] = unique(reshape(elem_sorted(:,idx),[],e_dim),'rows');
0082 elem2face = reshape(ia,[],e_dim+1);
0083 face2elem = calc_face2elem(elem2face);
0084 
0085 function face2elem = calc_face2elem(elem2face)
0086 % This is easier to understand but very slow
0087 %     n_face = max(elem2face(:));
0088 %     face2elem = zeros(n_face,2);
0089 %     for i= 1:n_face
0090 %         [el jnk] = find(elem2face==i);
0091 %         if numel(el)==1, el(2) = 0; end
0092 %         face2elem(i,:) = el;
0093 %     end
0094 %     bck = face2elem; face2elem=[];
0095     [n_elems, el_faces] = size(elem2face);
0096     elem2faceno = (1:n_elems)'*ones(1,el_faces);
0097     elem2faceno = elem2faceno(:);
0098     elem2face   = elem2face(:);
0099     face2elem(elem2face,2) = elem2faceno;
0100     % flipping will give us the other element for shared faces
0101     elem2faceno = flipud(elem2faceno);
0102     elem2face   = flipud(elem2face);
0103     face2elem(elem2face,1) = elem2faceno;
0104     % replace with zeros repeated entries (boundary faces)
0105     face2elem( face2elem(:,1) == face2elem(:,2), 2) = 0;
0106 
0107 % This function is obsolete
0108 function elem2face = calc_elem2face(face2elem, face_per_elem)
0109     n_elem = max(face2elem(:));
0110     elem2face = zeros(n_elem,face_per_elem);
0111     for i = 1:n_elem
0112         [f jnk] = find(face2elem==i);
0113         elem2face(i,:) = f;
0114     end
0115     
0116 function normals = calc_normals(mdl)
0117     [n_faces face_dim] = size(mdl.faces);
0118     switch face_dim
0119         case 2
0120             A = mdl.nodes(mdl.faces(:,1),:);
0121             B = mdl.nodes(mdl.faces(:,2),:);
0122             normals = (B-A)*[0 1; -1 0];
0123         case 3
0124             % vectorise cross product
0125             x1 = mdl.nodes(mdl.faces(:,2),1) - mdl.nodes(mdl.faces(:,1),1);
0126             y1 = mdl.nodes(mdl.faces(:,2),2) - mdl.nodes(mdl.faces(:,1),2);
0127             z1 = mdl.nodes(mdl.faces(:,2),3) - mdl.nodes(mdl.faces(:,1),3);
0128             x2 = mdl.nodes(mdl.faces(:,3),1) - mdl.nodes(mdl.faces(:,1),1);
0129             y2 = mdl.nodes(mdl.faces(:,3),2) - mdl.nodes(mdl.faces(:,1),2);
0130             z2 = mdl.nodes(mdl.faces(:,3),3) - mdl.nodes(mdl.faces(:,1),3);
0131             %(a2b3 ? a3b2, a3b1 ? a1b3, a1b2 ? a2b1).
0132             normals = zeros(n_faces,3);
0133             normals(:,1) = y1.*z2 - z1.*y2;
0134             normals(:,2) = z1.*x2 - x1.*z2;
0135             normals(:,3) = x1.*y2 - y1.*x2;
0136         otherwise;
0137             error('not 2D or 3D')
0138     end
0139     normals = normals./ repmat(sqrt(sum(normals.^2,2))',face_dim,1)';
0140     
0141     
0142 function len = calc_longest_edge(elems,nodes)
0143     [E_num E_dim] = size(elems);
0144 
0145     pairs = nchoosek(1:E_dim,2);
0146     len = zeros(E_num,1);
0147     for i = 1:size(pairs,1)
0148         a = nodes(elems(:,pairs(i,1)),:);
0149         b = nodes(elems(:,pairs(i,2)),:);
0150         tmp = sqrt(sum((a-b).^2,2));
0151         len = max(len,tmp);  
0152     end
0153     
0154 function do_unit_test
0155     % square
0156     nodes = [0 0; 0 1; 1 1; 1 0];
0157     elems = [1 2 3; 1 3 4];
0158     mdl = eidors_obj('fwd_model','square','nodes', nodes, 'elems', elems);
0159     out = fix_model(mdl);
0160     out.faces
0161     out.elem2face
0162     out.face2elem
0163 
0164     %cube
0165     nodes = [0 0 0; 0 1 0; 1 1 0; 1 0 0;...
0166              0 0 1; 0 1 1; 1 1 1; 1 0 1];
0167     elems = [1 2 3 6; 3 6 7 8; 1 5 6 8; 1 3 4 8; 1 3 6 8];
0168     mdl = eidors_obj('fwd_model','cube','nodes', nodes, 'elems', elems);
0169     out = fix_model(mdl);
0170     out.faces
0171     out.elem2face
0172     out.face2elem    
0173     
0174     mdl = mk_common_model('n3r2',16); mdl= mdl.fwd_model;
0175     out = fix_model(mdl);
0176