Description of fix

0001 function [mdl] = fix_model(mdl,opt)
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 - a struct with logical values specifying which fields to
0010 %        compute. Defaults to false for absent fields.
0011 %
0012 % Run fix_model('options') to get an all-false options struct, or
0013 % fix_model('options',true) to get an all-true one.
0014 % mdl.faces is only replaced if necessary. mdl.boundary is never replaced.
0015 %
0016 % OUTPUT:
0017 %    mdl - a copy of the input model with these additional fields:
0018 %       .boundary
0019 %       .boundary_face
0020 %       .faces
0021 %       .face2elem
0022 %       .elem2face
0023 %       .elem_centre
0024 %       .face_centre
0025 %       .face_area
0026 %       .normals
0027 %       .max_edge_len (per elem)
0028 %       .edges
0029 %       .edge_length
0030 %       .edge2elem
0031 %       .elem2edge
0032 %       .face2edge
0033 %       .node2elem
0034 %
0035 % For triangular meshes, edges and faces are the same.
0036 %
0037 % The elems will be reordered so that all faces are counter-clockwise.
0038 
0039 % (C) 2011 Bartlomiej Grychtol. Licensed under GPL v2 or v3
0040 % $Id: fix_model.m 5822 2018-08-27 13:03:01Z aadler $
0041 
0042 if ischar(mdl) && strcmp(mdl,'UNIT_TEST'); do_unit_test; return; end
0043 
0044 if ischar(mdl) && strcmp(mdl,'options'); 
0045    if nargin < 2, opt = false; end
0046    mdl = list_options(opt); 
0047    return 
0048 end
0049 doall = false;
0050 if nargin > 1
0051    opt = fix_options(mdl,opt);
0052 else
0053    doall = true; opt=struct;
0054 end
0055 
0056 % prepare a model with geometry only
0057 tmp.nodes = mdl.nodes;
0058 tmp.elems = mdl.elems;
0059 tmp.type  = mdl.type;
0060 try, tmp.boundary = mdl.boundary; end
0061 try, tmp.faces = mdl.faces; end
0062 
0063 copt.cache_obj = {tmp, doall, opt};
0064 copt.fstr      = 'fix_model';
0065 
0066 tmp = eidors_cache( @do_fix_model, {tmp, doall, opt}, copt);
0067 
0068 % copy new fields to mdl
0069 flds = fieldnames(tmp); flds(1:3) = []; % ignore nodes, elems and type
0070 for i = 1:numel(flds)
0071    mdl.(flds{i}) = tmp.(flds{i});
0072 end
0073 
0074 % This code may not stay here. Think about where to put it
0075 function fmdl = remove_unused_nodes( fmdl );
0076    usednodes = unique(fmdl.elems(:));
0077    if max(usednodes) > num_nodes(fmdl)
0078       error('more nodes are used than exist');
0079    end
0080    nidx = zeros(num_nodes(fmdl),1);
0081    nidx(usednodes) = 1:length(usednodes);
0082    fmdl.nodes(nidx==0,:) = [];
0083    fmdl.elems = reshape(nidx(fmdl.elems),[],size(fmdl.elems,2));
0084 
0085    for i=1:length(fmdl.electrode)
0086       fmdl.electrode(i).nodes =  nidx( ...
0087          fmdl.electrode(i).nodes);
0088    end
0089    fmdl.boundary = find_boundary(fmdl);
0090    fmdl.gnd_node = usednodes(fmdl.gnd_node);
0091 
0092 
0093 % Complete the function
0094 function mdl = do_fix_model(mdl, doall, opt)
0095 
0096 if doall || opt.boundary
0097    if ~isfield(mdl,'boundary')
0098       mdl.boundary = find_boundary(mdl);
0099    end
0100 end
0101 if doall || opt.faces
0102    if elem_dim(mdl) == mdl_dim(mdl);
0103       [mdl.faces mdl.elem2face] = calc_faces(mdl);
0104    else
0105       % surface mesh
0106       mdl.faces = sort(mdl.elems,2);
0107       mdl.elem2face = (1:length(mdl.faces))';
0108    end
0109 end
0110 if doall || opt.face2elem
0111     mdl.face2elem = calc_face2elem(mdl.elem2face);
0112 end
0113 if doall || opt.boundary_face
0114    mdl.boundary_face = mdl.face2elem(:,2)==0;
0115 end
0116 if doall || opt.elem_centre
0117    mdl.elem_centre = interp_mesh(mdl, 0);
0118 end
0119 if doall || opt.face_centre
0120    tmp = mdl;
0121    tmp.elems = tmp.faces;
0122    mdl.face_centre = interp_mesh(tmp,0);
0123 end
0124 if doall || opt.max_edge_len
0125    mdl.max_edge_len = calc_longest_edge(mdl.elems,mdl.nodes);
0126 end
0127 if doall || opt.elem_volume
0128    mdl.elem_volume = get_elem_volume(mdl);
0129 end
0130 if doall || opt.face_area
0131    if elem_dim(mdl) == 2
0132       mdl.face_area = get_elem_volume(mdl);
0133    else
0134       mdl.face_area = calc_face_area(mdl);
0135    end
0136 end
0137 el_dim = elem_dim(mdl);
0138 if doall || opt.edges
0139     if mdl_dim(mdl)==3 %that's unlikely to work for higher order elements
0140         [mdl.edges mdl.elem2edge] = calc_faces(mdl,2);
0141     else 
0142         mdl.edges = mdl.faces;
0143         mdl.elem2edge = mdl.elem2face;
0144     end
0145 end
0146 if doall || opt.node2elem
0147     mdl.node2elem = calc_edge2elem(mdl.elems,size(mdl.nodes,1));
0148 end
0149 if doall || opt.edge2elem
0150     mdl.edge2elem = calc_edge2elem(mdl.elem2edge, size(mdl.edges,1));
0151 end
0152 
0153 if doall || opt.face2edge
0154    if el_dim <3
0155       mdl.face2edge = mdl.elem2edge;
0156    else
0157       mdl.face2edge = uint32(calc_face2edge(mdl));
0158    end
0159 end
0160 
0161 if doall || opt.edge_length
0162    mdl.edge_length = calc_edge_length(mdl);
0163 end
0164 
0165 if doall || opt.linear_reorder
0166    mdl = linear_reorder(mdl); %counter-clockwise
0167 end
0168 if doall || opt.normals
0169    mdl.normals = calc_normals(mdl);
0170 end
0171 if doall || opt.inner_normal
0172    if mdl_dim(mdl) == elem_dim(mdl);
0173       mdl.inner_normal = test_inner_normal( mdl );
0174    else
0175       eidors_msg('@@@ Inner normal test for surface meshes not implemented.',1);
0176    end
0177 end
0178 
0179 % decrease memory footprint
0180 mdl.elems = uint32(mdl.elems);
0181 if doall || opt.faces
0182    mdl.faces = uint32(mdl.faces);
0183    mdl.elem2face = uint32(mdl.elem2face);
0184 end
0185 if doall || opt.face2elem
0186    mdl.face2elem = uint32(mdl.face2elem);
0187 end     
0188 if doall || opt.edges
0189    mdl.edges = uint32(mdl.edges);
0190    mdl.elem2edge = uint32(mdl.elem2edge);
0191 end
0192 if doall || opt.edge2elem
0193    mdl.edge2elem = logical(mdl.edge2elem);
0194 end     
0195 
0196 % Test whether normal points into or outside
0197 % mdl.inner_normal(i,j) = 1 if face j of elem i points in
0198 function inner_normal = test_inner_normal( mdl )
0199    inner_normal = false(size(mdl.elem2face));
0200    d = elem_dim(mdl) + 1;
0201    for i=1:num_elems(mdl);
0202       el_faces = mdl.elem2face(i,:);
0203       el_ctr   = repmat( mdl.elem_centre(i,:), d, 1);
0204       vec_fa_el= el_ctr -  mdl.face_centre(el_faces,:);
0205       normal_i  = mdl.normals(el_faces,:);
0206       dot_prod = sum( normal_i.*vec_fa_el, 2 );
0207       inner_normal(i,:) = dot_prod' > 0;
0208    end
0209 
0210 function [faces,  elem2face] = calc_faces(mdl, facedim)
0211 
0212 e_dim = elem_dim(mdl);
0213 if nargin == 1
0214     facedim = e_dim;
0215 end
0216 
0217 idx = nchoosek(1:e_dim+1, facedim);
0218 elem_sorted = sort(mdl.elems,2);
0219 [faces ib ia] = unique(reshape(elem_sorted(:,idx),[],facedim),'rows');
0220 elem2face = reshape(ia,[],size(idx,1));
0221 
0222 function edge2elem = calc_edge2elem(elem2edge,n_edges)
0223 
0224     [n_elems, el_faces] = size(elem2edge);
0225     elem2edgeno = (1:n_elems)'*ones(1,el_faces);
0226     elem2edgeno = elem2edgeno(:);
0227     elem2edge   = elem2edge(:);
0228     edge2elem = sparse(elem2edge,elem2edgeno,ones(size(elem2edgeno)),n_edges,n_elems);
0229 
0230 function f2e = calc_face2edge(mdl)
0231 %faces and edges are both row-wise sorted
0232 nf = length(mdl.faces);
0233 list(1:3:3*nf,:) = mdl.faces(:,1:2);
0234 list(2:3:3*nf,:) = mdl.faces(:,[1 3]);
0235 list(3:3:3*nf,:) = mdl.faces(:,2:3);
0236 [jnk f2e] = ismember(list, mdl.edges, 'rows');
0237 f2e = reshape(f2e,3,[])';
0238 
0239 function face2elem = calc_face2elem(elem2face)
0240 % This is easier to understand but very slow
0241 %     n_face = max(elem2face(:));
0242 %     face2elem = zeros(n_face,2);
0243 %     for i= 1:n_face
0244 %         [el jnk] = find(elem2face==i);
0245 %         if numel(el)==1, el(2) = 0; end
0246 %         face2elem(i,:) = el;
0247 %     end
0248 %     bck = face2elem; face2elem=[];
0249     [n_elems, el_faces] = size(elem2face);
0250     elem2faceno = (1:n_elems)'*ones(1,el_faces);
0251     elem2faceno = elem2faceno(:);
0252     elem2face   = elem2face(:);
0253     face2elem(elem2face,2) = elem2faceno;
0254     % flipping will give us the other element for shared faces
0255     elem2faceno = flipud(elem2faceno);
0256     elem2face   = flipud(elem2face);
0257     face2elem(elem2face,1) = elem2faceno;
0258     % replace with zeros repeated entries (boundary faces)
0259     face2elem( face2elem(:,1) == face2elem(:,2), 2) = 0;
0260 
0261 % This function is obsolete
0262 function elem2face = calc_elem2face(face2elem, face_per_elem)
0263     n_elem = max(face2elem(:));
0264     elem2face = zeros(n_elem,face_per_elem);
0265     for i = 1:n_elem
0266         [f jnk] = find(face2elem==i);
0267         elem2face(i,:) = f;
0268     end
0269     
0270 function normals = calc_normals(mdl)
0271     [n_faces face_dim] = size(mdl.faces);
0272     switch face_dim
0273         case 2
0274             A = mdl.nodes(mdl.faces(:,1),:);
0275             B = mdl.nodes(mdl.faces(:,2),:);
0276             normals = (B-A)*[0 1; -1 0];
0277         case 3
0278             % vectorise cross product
0279             x1 = mdl.nodes(mdl.faces(:,2),1) - mdl.nodes(mdl.faces(:,1),1);
0280             y1 = mdl.nodes(mdl.faces(:,2),2) - mdl.nodes(mdl.faces(:,1),2);
0281             z1 = mdl.nodes(mdl.faces(:,2),3) - mdl.nodes(mdl.faces(:,1),3);
0282             x2 = mdl.nodes(mdl.faces(:,3),1) - mdl.nodes(mdl.faces(:,1),1);
0283             y2 = mdl.nodes(mdl.faces(:,3),2) - mdl.nodes(mdl.faces(:,1),2);
0284             z2 = mdl.nodes(mdl.faces(:,3),3) - mdl.nodes(mdl.faces(:,1),3);
0285             %(a2b3 ? a3b2, a3b1 ? a1b3, a1b2 ? a2b1).
0286             normals = zeros(n_faces,3);
0287             normals(:,1) = y1.*z2 - z1.*y2;
0288             normals(:,2) = z1.*x2 - x1.*z2;
0289             normals(:,3) = x1.*y2 - y1.*x2;
0290         otherwise;
0291             error('not 2D or 3D')
0292     end
0293     normals = normals./ repmat(sqrt(sum(normals.^2,2))',face_dim,1)';
0294     
0295  function A = calc_face_area(mdl)
0296 A = mdl.nodes(mdl.faces(:,2),:) - mdl.nodes(mdl.faces(:,1),:);
0297 B = mdl.nodes(mdl.faces(:,3),:) - mdl.nodes(mdl.faces(:,1),:);
0298 A = sqrt(sum(cross3(A,B).^2,2))/2;
0299 
0300 function L = calc_edge_length(mdl)
0301 L = sqrt(sum( (mdl.nodes(mdl.edges(:,1),:) ...
0302              - mdl.nodes(mdl.edges(:,2),:) ).^2 ,2 ));
0303     
0304 function len = calc_longest_edge(elems,nodes)
0305     [E_num E_dim] = size(elems);
0306 
0307     pairs = nchoosek(1:E_dim,2);
0308     len = zeros(E_num,1);
0309     for i = 1:size(pairs,1)
0310         a = nodes(elems(:,pairs(i,1)),:);
0311         b = nodes(elems(:,pairs(i,2)),:);
0312         tmp = sqrt(sum((a-b).^2,2));
0313         len = max(len,tmp);  
0314     end
0315     
0316 function out = fix_options(mdl, opt)
0317     out = list_options(false);
0318 %     out.linear_reorder = 1;
0319     flds = fieldnames(opt);
0320     for i = 1:length(flds)
0321        try
0322        out.(flds{i}) = opt.(flds{i});
0323        catch
0324           warning(sprintf('Option %s not recognised. Ignoring', flds{i}));
0325        end
0326     end 
0327     if out.boundary_face
0328        out.face2elem = true;
0329     end
0330     if out.inner_normal
0331        out.normals = true;
0332        out.elem_centre = true;
0333        out.face_centre = true;
0334     end
0335     if any([ out.boundary_face out.face_centre out.normals]) && ~isfield(mdl,'faces')
0336           out.faces = true;
0337     end
0338     if any([out.face2elem out.elem2face out.face_area])
0339        out.faces = true;
0340     end
0341     if any([out.edge2elem out.elem2edge out.edge_length out.face2edge])
0342         out.edges = true;
0343     end
0344     if out.edges && elem_dim(mdl) < 4
0345         out.faces = true;
0346     end
0347 
0348     
0349 function out = list_options(val)
0350     nodes = [0 0; 0 1; 1 1; 1 0];
0351     elems = [1 2 3; 1 3 4];
0352     mdl = eidors_obj('fwd_model','square','nodes', nodes, 'elems', elems);
0353     out = fix_model(mdl);
0354     out = rmfield(out,{'elems','nodes','name','type'});
0355     flds = fieldnames(out);
0356     for i = 1:length(flds)
0357        out.(flds{i}) = val;
0358     end
0359     out.linear_reorder = 0;
0360     
0361 function do_unit_test
0362     % square
0363     nodes = [0 0; 0 1; 1 1; 1 0];
0364     elems = [1 2 3; 1 3 4];
0365     mdl = eidors_obj('fwd_model','square','nodes', nodes, 'elems', elems);
0366     out = fix_model(mdl);
0367     unit_test_cmp('2d: faces'    ,out.faces    ,[1,2;1,3;1,4;2,3;3,4]);
0368     unit_test_cmp('2d: elem2face',out.elem2face,[1,2,4;2,3,5]);
0369     unit_test_cmp('2d: face2elem',out.face2elem,[1,0;2,1;2,0;1,0;2,0]);
0370 
0371     %cube
0372     nodes = [0 0 0; 0 1 0; 1 1 0; 1 0 0;...
0373              0 0 1; 0 1 1; 1 1 1; 1 0 1];
0374     elems = [1 2 3 6; 3 6 7 8; 1 5 6 8; 1 3 4 8; 1 3 6 8];
0375     mdl = eidors_obj('fwd_model','cube','nodes', nodes, 'elems', elems);
0376     out = fix_model(mdl);
0377     unit_test_cmp('3d: faces'    ,out.faces(1:4,:), [1,2,3;1,2,6;1,3,4;1,3,6]);
0378     unit_test_cmp('3d: elem2face',out.elem2face, [1,2,4,10; 12,13,14,16;
0379             7,8,9,15; 3,5,6,11; 4,5,9,13]);
0380     unit_test_cmp('3d: face2elem',out.face2elem(1:5,:), [1,0; 1,0; 4,0; 5,1; 4,5]);
0381 
0382     
0383     % test options
0384     opt = fix_model('options',false);
0385     flds = fieldnames(opt);
0386     % if there are no errors, option interdependence is dealt with
0387     % correctly
0388     for i = 1:length(flds)
0389        opt.(flds{i}) = true;
0390        out = fix_model(mdl, opt);
0391        opt.(flds{i}) = false;
0392     end
0393     
0394     mdl = mk_common_model('n3r2',[16,2]); mdl= mdl.fwd_model;
0395     out = fix_model(mdl);
0396
fix_model

PURPOSE

SYNOPSIS

DESCRIPTION

CROSS-REFERENCE INFORMATION

SUBFUNCTIONS

SOURCE CODE
