0001 function [mdl] = fix_model(mdl,opt)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
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
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
0069 flds = fieldnames(tmp); flds(1:3) = [];
0070 for i = 1:numel(flds)
0071 mdl.(flds{i}) = tmp.(flds{i});
0072 end
0073
0074
0075
0076
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
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
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);
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
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
0180
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
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
0224
0225
0226
0227
0228
0229
0230
0231
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
0238 elem2faceno = flipud(elem2faceno);
0239 elem2face = flipud(elem2face);
0240 face2elem(elem2face,1) = elem2faceno;
0241
0242 face2elem( face2elem(:,1) == face2elem(:,2), 2) = 0;
0243
0244
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
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
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
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
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
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
0367 opt = fix_model('options',false);
0368 flds = fieldnames(opt);
0369
0370
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