fix_model

PURPOSE ^

FIX_MODEL: Add useful fields to a model

SYNOPSIS ^

function [mdl] = fix_model(mdl,opt)

DESCRIPTION ^

 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 - a struct with logical values specifying which fields to
        compute. Defaults to false for absent fields. 

 Run fix_model('options') to get an all-false options struct, or
 fix_model('options',true) to get an all-true one.
 mdl.faces is only replaced if necessary. mdl.boundary is never replaced.

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

 For triangular meshes, edges and faces are the same.

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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 5361 2017-03-16 12:40:08Z bgrychtol $
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 
0075 
0076 % Complete the function
0077 function mdl = do_fix_model(mdl, doall, opt)
0078 
0079 if doall || opt.boundary
0080    if ~isfield(mdl,'boundary')
0081       mdl.boundary = find_boundary(mdl);
0082    end
0083 end
0084 if doall || opt.faces
0085    if elem_dim(mdl) == mdl_dim(mdl);
0086       [mdl.faces mdl.elem2face] = calc_faces(mdl);
0087    else
0088       % surface mesh
0089       mdl.faces = sort(mdl.elems,2);
0090       mdl.elem2face = (1:length(mdl.faces))';
0091    end
0092 end
0093 if doall || opt.face2elem
0094     mdl.face2elem = calc_face2elem(mdl.elem2face);
0095 end
0096 if doall || opt.boundary_face
0097    mdl.boundary_face = mdl.face2elem(:,2)==0;
0098 end
0099 if doall || opt.elem_centre
0100    mdl.elem_centre = interp_mesh(mdl, 0);
0101 end
0102 if doall || opt.face_centre
0103    tmp = mdl;
0104    tmp.elems = tmp.faces;
0105    mdl.face_centre = interp_mesh(tmp,0);
0106 end
0107 if doall || opt.max_edge_len
0108    mdl.max_edge_len = calc_longest_edge(mdl.elems,mdl.nodes);
0109 end
0110 if doall || opt.elem_volume
0111    mdl.elem_volume = get_elem_volume(mdl);
0112 end
0113 if doall || opt.face_area
0114    if elem_dim(mdl) == 2
0115       mdl.face_area = get_elem_volume(mdl);
0116    else
0117       mdl.face_area = calc_face_area(mdl);
0118    end
0119 end
0120 el_dim = elem_dim(mdl);
0121 if doall || opt.edges
0122     if mdl_dim(mdl)==3 %that's unlikely to work for higher order elements
0123         [mdl.edges mdl.elem2edge] = calc_faces(mdl,2);
0124     else 
0125         mdl.edges = mdl.faces;
0126         mdl.elem2edge = mdl.elem2face;
0127     end
0128 end
0129 if doall || opt.node2elem
0130     mdl.node2elem = calc_edge2elem(mdl.elems,size(mdl.nodes,1));
0131 end
0132 if doall || opt.edge2elem
0133     mdl.edge2elem = calc_edge2elem(mdl.elem2edge, size(mdl.edges,1));
0134 end
0135 
0136 if doall || opt.face2edge
0137    if el_dim <3
0138       mdl.face2edge = mdl.elem2edge;
0139    else
0140       mdl.face2edge = uint32(calc_face2edge(mdl));
0141    end
0142 end
0143 
0144 if doall || opt.edge_length
0145    mdl.edge_length = calc_edge_length(mdl);
0146 end
0147 
0148 if doall || opt.linear_reorder
0149    mdl = linear_reorder(mdl); %counter-clockwise
0150 end
0151 if doall || opt.normals
0152    mdl.normals = calc_normals(mdl);
0153 end
0154 if doall || opt.inner_normal
0155    if mdl_dim(mdl) == elem_dim(mdl);
0156       mdl.inner_normal = test_inner_normal( mdl );
0157    else
0158       eidors_msg('@@@ Inner normal test for surface meshes not implemented.',1);
0159    end
0160 end
0161 
0162 % decrease memory footprint
0163 mdl.elems = uint32(mdl.elems);
0164 if doall || opt.faces
0165    mdl.faces = uint32(mdl.faces);
0166    mdl.elem2face = uint32(mdl.elem2face);
0167 end
0168 if doall || opt.face2elem
0169    mdl.face2elem = uint32(mdl.face2elem);
0170 end     
0171 if doall || opt.edges
0172    mdl.edges = uint32(mdl.edges);
0173    mdl.elem2edge = uint32(mdl.elem2edge);
0174 end
0175 if doall || opt.edge2elem
0176    mdl.edge2elem = logical(mdl.edge2elem);
0177 end     
0178 
0179 % Test whether normal points into or outside
0180 % mdl.inner_normal(i,j) = 1 if face j of elem i points in
0181 function inner_normal = test_inner_normal( mdl )
0182    inner_normal = false(size(mdl.elem2face));
0183    d = elem_dim(mdl) + 1;
0184    for i=1:num_elems(mdl);
0185       el_faces = mdl.elem2face(i,:);
0186       el_ctr   = repmat( mdl.elem_centre(i,:), d, 1);
0187       vec_fa_el= el_ctr -  mdl.face_centre(el_faces,:);
0188       normal_i  = mdl.normals(el_faces,:);
0189       dot_prod = sum( normal_i.*vec_fa_el, 2 );
0190       inner_normal(i,:) = dot_prod' > 0;
0191    end
0192 
0193 function [faces,  elem2face] = calc_faces(mdl, facedim)
0194 
0195 e_dim = elem_dim(mdl);
0196 if nargin == 1
0197     facedim = e_dim;
0198 end
0199 
0200 idx = nchoosek(1:e_dim+1, facedim);
0201 elem_sorted = sort(mdl.elems,2);
0202 [faces ib ia] = unique(reshape(elem_sorted(:,idx),[],facedim),'rows');
0203 elem2face = reshape(ia,[],size(idx,1));
0204 
0205 function edge2elem = calc_edge2elem(elem2edge,n_edges)
0206 
0207     [n_elems, el_faces] = size(elem2edge);
0208     elem2edgeno = (1:n_elems)'*ones(1,el_faces);
0209     elem2edgeno = elem2edgeno(:);
0210     elem2edge   = elem2edge(:);
0211     edge2elem = sparse(elem2edge,elem2edgeno,ones(size(elem2edgeno)),n_edges,n_elems);
0212 
0213 function f2e = calc_face2edge(mdl)
0214 %faces and edges are both row-wise sorted
0215 nf = length(mdl.faces);
0216 list(1:3:3*nf,:) = mdl.faces(:,1:2);
0217 list(2:3:3*nf,:) = mdl.faces(:,[1 3]);
0218 list(3:3:3*nf,:) = mdl.faces(:,2:3);
0219 [jnk f2e] = ismember(list, mdl.edges, 'rows');
0220 f2e = reshape(f2e,3,[])';
0221 
0222 function face2elem = calc_face2elem(elem2face)
0223 % This is easier to understand but very slow
0224 %     n_face = max(elem2face(:));
0225 %     face2elem = zeros(n_face,2);
0226 %     for i= 1:n_face
0227 %         [el jnk] = find(elem2face==i);
0228 %         if numel(el)==1, el(2) = 0; end
0229 %         face2elem(i,:) = el;
0230 %     end
0231 %     bck = face2elem; face2elem=[];
0232     [n_elems, el_faces] = size(elem2face);
0233     elem2faceno = (1:n_elems)'*ones(1,el_faces);
0234     elem2faceno = elem2faceno(:);
0235     elem2face   = elem2face(:);
0236     face2elem(elem2face,2) = elem2faceno;
0237     % flipping will give us the other element for shared faces
0238     elem2faceno = flipud(elem2faceno);
0239     elem2face   = flipud(elem2face);
0240     face2elem(elem2face,1) = elem2faceno;
0241     % replace with zeros repeated entries (boundary faces)
0242     face2elem( face2elem(:,1) == face2elem(:,2), 2) = 0;
0243 
0244 % This function is obsolete
0245 function elem2face = calc_elem2face(face2elem, face_per_elem)
0246     n_elem = max(face2elem(:));
0247     elem2face = zeros(n_elem,face_per_elem);
0248     for i = 1:n_elem
0249         [f jnk] = find(face2elem==i);
0250         elem2face(i,:) = f;
0251     end
0252     
0253 function normals = calc_normals(mdl)
0254     [n_faces face_dim] = size(mdl.faces);
0255     switch face_dim
0256         case 2
0257             A = mdl.nodes(mdl.faces(:,1),:);
0258             B = mdl.nodes(mdl.faces(:,2),:);
0259             normals = (B-A)*[0 1; -1 0];
0260         case 3
0261             % vectorise cross product
0262             x1 = mdl.nodes(mdl.faces(:,2),1) - mdl.nodes(mdl.faces(:,1),1);
0263             y1 = mdl.nodes(mdl.faces(:,2),2) - mdl.nodes(mdl.faces(:,1),2);
0264             z1 = mdl.nodes(mdl.faces(:,2),3) - mdl.nodes(mdl.faces(:,1),3);
0265             x2 = mdl.nodes(mdl.faces(:,3),1) - mdl.nodes(mdl.faces(:,1),1);
0266             y2 = mdl.nodes(mdl.faces(:,3),2) - mdl.nodes(mdl.faces(:,1),2);
0267             z2 = mdl.nodes(mdl.faces(:,3),3) - mdl.nodes(mdl.faces(:,1),3);
0268             %(a2b3 ? a3b2, a3b1 ? a1b3, a1b2 ? a2b1).
0269             normals = zeros(n_faces,3);
0270             normals(:,1) = y1.*z2 - z1.*y2;
0271             normals(:,2) = z1.*x2 - x1.*z2;
0272             normals(:,3) = x1.*y2 - y1.*x2;
0273         otherwise;
0274             error('not 2D or 3D')
0275     end
0276     normals = normals./ repmat(sqrt(sum(normals.^2,2))',face_dim,1)';
0277     
0278  function A = calc_face_area(mdl)
0279 A = mdl.nodes(mdl.faces(:,2),:) - mdl.nodes(mdl.faces(:,1),:);
0280 B = mdl.nodes(mdl.faces(:,3),:) - mdl.nodes(mdl.faces(:,1),:);
0281 A = sqrt(sum(cross3(A,B).^2,2))/2;
0282 
0283 function L = calc_edge_length(mdl)
0284 L = sqrt(sum( (mdl.nodes(mdl.edges(:,1),:) ...
0285              - mdl.nodes(mdl.edges(:,2),:) ).^2 ,2 ));
0286     
0287 function len = calc_longest_edge(elems,nodes)
0288     [E_num E_dim] = size(elems);
0289 
0290     pairs = nchoosek(1:E_dim,2);
0291     len = zeros(E_num,1);
0292     for i = 1:size(pairs,1)
0293         a = nodes(elems(:,pairs(i,1)),:);
0294         b = nodes(elems(:,pairs(i,2)),:);
0295         tmp = sqrt(sum((a-b).^2,2));
0296         len = max(len,tmp);  
0297     end
0298     
0299 function out = fix_options(mdl, opt)
0300     out = list_options(false);
0301 %     out.linear_reorder = 1;
0302     flds = fieldnames(opt);
0303     for i = 1:length(flds)
0304        try
0305        out.(flds{i}) = opt.(flds{i});
0306        catch
0307           warning(sprintf('Option %s not recognised. Ignoring', flds{i}));
0308        end
0309     end 
0310     if out.boundary_face
0311        out.face2elem = true;
0312     end
0313     if out.inner_normal
0314        out.normals = true;
0315        out.elem_centre = true;
0316        out.face_centre = true;
0317     end
0318     if any([ out.boundary_face out.face_centre out.normals]) && ~isfield(mdl,'faces')
0319           out.faces = true;
0320     end
0321     if any([out.face2elem out.elem2face out.face_area])
0322        out.faces = true;
0323     end
0324     if any([out.edge2elem out.elem2edge out.edge_length out.face2edge])
0325         out.edges = true;
0326     end
0327     if out.edges && elem_dim(mdl) < 4
0328         out.faces = true;
0329     end
0330 
0331     
0332 function out = list_options(val)
0333     nodes = [0 0; 0 1; 1 1; 1 0];
0334     elems = [1 2 3; 1 3 4];
0335     mdl = eidors_obj('fwd_model','square','nodes', nodes, 'elems', elems);
0336     out = fix_model(mdl);
0337     out = rmfield(out,{'elems','nodes','name','type'});
0338     flds = fieldnames(out);
0339     for i = 1:length(flds)
0340        out.(flds{i}) = val;
0341     end
0342     out.linear_reorder = 0;
0343     
0344 function do_unit_test
0345     % square
0346     nodes = [0 0; 0 1; 1 1; 1 0];
0347     elems = [1 2 3; 1 3 4];
0348     mdl = eidors_obj('fwd_model','square','nodes', nodes, 'elems', elems);
0349     out = fix_model(mdl);
0350     unit_test_cmp('2d: faces'    ,out.faces    ,[1,2;1,3;1,4;2,3;3,4]);
0351     unit_test_cmp('2d: elem2face',out.elem2face,[1,2,4;2,3,5]);
0352     unit_test_cmp('2d: face2elem',out.face2elem,[1,0;2,1;2,0;1,0;2,0]);
0353 
0354     %cube
0355     nodes = [0 0 0; 0 1 0; 1 1 0; 1 0 0;...
0356              0 0 1; 0 1 1; 1 1 1; 1 0 1];
0357     elems = [1 2 3 6; 3 6 7 8; 1 5 6 8; 1 3 4 8; 1 3 6 8];
0358     mdl = eidors_obj('fwd_model','cube','nodes', nodes, 'elems', elems);
0359     out = fix_model(mdl);
0360     unit_test_cmp('3d: faces'    ,out.faces(1:4,:), [1,2,3;1,2,6;1,3,4;1,3,6]);
0361     unit_test_cmp('3d: elem2face',out.elem2face, [1,2,4,10; 12,13,14,16;
0362             7,8,9,15; 3,5,6,11; 4,5,9,13]);
0363     unit_test_cmp('3d: face2elem',out.face2elem(1:5,:), [1,0; 1,0; 4,0; 5,1; 4,5]);
0364 
0365     
0366     % test options
0367     opt = fix_model('options',false);
0368     flds = fieldnames(opt);
0369     % if there are no errors, option interdependence is dealt with
0370     % correctly
0371     for i = 1:length(flds)
0372        opt.(flds{i}) = true;
0373        out = fix_model(mdl, opt);
0374        opt.(flds{i}) = false;
0375     end
0376     
0377     mdl = mk_common_model('n3r2',[16,2]); mdl= mdl.fwd_model;
0378     out = fix_model(mdl);
0379

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005