0001 function [Q mdl] = calc_mesh_quality(mdl, show)
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
0043
0044
0045
0046
0047
0048
0049
0050
0051
0052
0053
0054
0055
0056
0057
0058
0059
0060
0061
0062
0063
0064
0065
0066
0067
0068
0069
0070
0071
0072
0073
0074
0075
0076
0077
0078
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105 if ischar(mdl) && strcmp(mdl, 'UNIT_TEST'), do_unit_test; return; end
0106
0107 if nargin < 2
0108 show = 0;
0109 end
0110
0111 opt.cache_obj = {mdl.nodes, mdl.elems};
0112 opt.fstr = 'calc_mesh_quality';
0113 [Q mdl] = eidors_cache(@do_quality_calc,{mdl},opt);
0114
0115 if show
0116 display_figs(Q)
0117 end
0118
0119
0120 function [Q mdl] = do_quality_calc(mdl)
0121
0122 opt.elem2edge = 1;
0123 opt.face2elem = 1;
0124 opt.boundary = 1;
0125 opt.boundary_face = 1;
0126 opt.normals = 1;
0127 opt.face_centre = 1;
0128 opt.inner_normal = 1;
0129 opt.face2edge = 1;
0130 opt.face_area = 1;
0131 opt.edge_length = 1;
0132 mdl = fix_model(mdl, opt);
0133
0134
0135
0136 if elem_dim(mdl) == 3
0137 mdl.max_edge_length = max(mdl.edge_length(mdl.elem2edge),[],2);
0138 mdl.min_edge_length = min(mdl.edge_length(mdl.elem2edge),[],2);
0139 mdl.elem_vol = elem_volume(mdl);
0140 mdl.elem_area = elem_area(mdl);
0141 mdl.rad_insphere = insphere_radius(mdl);
0142 mdl.rad_circumsphere = circumsphere_radius(mdl);
0143 mdl.tet_alt = tet_altitudes(mdl);
0144
0145
0146
0147 mdl.solid_angle = solid_angles(mdl);
0148 mdl.dihedral_angle = dihedral_angles(mdl);
0149 end
0150
0151
0152 mdl.rad_incircle = incircle_radius(mdl);
0153 mdl.rad_circumcircle = circumcircle_radius(mdl);
0154 mdl.tri_alt = tri_altitudes(mdl);
0155 mdl.tri_angle = tri_angles(mdl);
0156
0157
0158 Q.tri.NSR = tri_normalized_shape_ratio(mdl);
0159 Q.tri.mu = tri_mu(mdl);
0160 Q.tri.eta = tri_eta(mdl);
0161 Q.tri.theta = tri_theta(mdl);
0162 Q.tri.iota = tri_iota(mdl);
0163 Q.tri.kappa = tri_kappa(mdl);
0164 Q.tri.min_angle = tri_min_angle(mdl);
0165
0166
0167 f = fieldnames(Q.tri);
0168 for i = 1:length(f)
0169 Q.tri.(f{i}) = Q.tri.(f{i})(mdl.boundary_face);
0170 end
0171
0172
0173 if elem_dim(mdl) == 3
0174
0175 Q.tet.NSR = tet_normalized_shape_ratio(mdl);
0176
0177
0178 Q.tet.mu = mu(mdl);
0179 Q.tet.tau = tau(mdl);
0180 Q.tet.reg = R_measure(mdl);
0181
0182
0183
0184 Q.tet.zeta = zeta(mdl);
0185 Q.tet.eta = eta(mdl);
0186 Q.tet.alpha = alpha(mdl);
0187 Q.tet.gamma = gamma(mdl);
0188 Q.tet.min_angle = tet_min_angle(mdl);
0189 end
0190
0191 function display_figs(Q);
0192 f = figure; set(f,'Name','Surface triangle quality');
0193 f = fieldnames(Q.tri);
0194 for i = 1:length(f)
0195 subplot(3,3,i)
0196 hist(Q.tri.(f{i}),100);
0197 xlabel(strrep(f{i},'_','\_'));
0198 axis tight
0199 xlim([0 1]);
0200 end
0201 if ~isfield(Q,'tet'), return, end;
0202 f = figure; set(f,'Name','Tetrahedron quality');
0203 f = fieldnames(Q.tet);
0204 for i = 1:length(f)
0205 subplot(3,3,i)
0206 hist(Q.tet.(f{i}),100);
0207 xlabel(strrep(f{i},'_','\_'));
0208 axis tight
0209 xlim([0 1]);
0210 end
0211
0212
0213
0214
0215
0216
0217 function a = tet_min_angle(mdl)
0218 a = min(mdl.dihedral_angle,[],2) / acos(1/3);
0219
0220 function a = tri_min_angle(mdl)
0221 a = 3*min(mdl.tri_angle,[],2)/pi;
0222
0223 function k = tri_kappa(mdl)
0224 k = 2/sqrt(3) * min(mdl.tri_alt,[],2) ...
0225 ./ max(mdl.edge_length(mdl.face2edge),[],2);
0226
0227 function i = tri_iota(mdl)
0228 i = 12*sqrt(3)*mdl.face_area ./ sum(mdl.edge_length(mdl.face2edge),2).^2;
0229
0230 function G = gamma(mdl)
0231 G = 6 * sqrt(2) * mdl.elem_vol ./ sqrt(mean(mdl.edge_length(mdl.elem2edge).^2,2)).^3 ;
0232
0233 function A = alpha(mdl)
0234 A = 6 * sqrt(2) * mdl.elem_vol ./ mean(mdl.edge_length(mdl.elem2edge),2).^3 ;
0235
0236 function E = eta(mdl)
0237 E = 12 * (3*mdl.elem_vol).^(2/3) ./ sum(mdl.edge_length(mdl.elem2edge).^2,2);
0238
0239 function E = tri_eta(mdl)
0240 E = 4 * sqrt(3) * mdl.face_area ./ sum(mdl.edge_length(mdl.face2edge).^2,2);
0241
0242 function Z = zeta(mdl)
0243 Z = 3^7 * mdl.elem_vol.^4 ./ ...
0244 sum(mdl.face_area(mdl.elem2face).^2 ,2).^3;
0245
0246 function S = S_measure(mdl)
0247 S = mdl.elem_vol ./ sum( mdl.edge_length(mdl.elem2edge).^2,2);
0248
0249 function R = R_star(mdl)
0250 R = sqrt(2)*mdl.elem_vol ./ sum(mdl.edge_length(mdl.elem2edge),2);
0251
0252 function R = R_measure(mdl)
0253 R = 4*mdl.rad_insphere ./ max(mdl.tet_alt,[],2);
0254
0255 function t = tau(mdl)
0256 t = mdl.min_edge_length ./ mdl.max_edge_length;
0257
0258 function m = mu(mdl)
0259 m = 2*sqrt(6)*mdl.rad_insphere./mdl.max_edge_length;
0260
0261 function m = tri_mu(mdl)
0262 m = 2*sqrt(3)*mdl.rad_incircle./ max(mdl.edge_length(mdl.face2edge),[],2);
0263
0264 function t = tri_theta(mdl)
0265 t = 2*min(mdl.tri_alt,[],2) ./ sqrt(sum(mdl.edge_length(mdl.face2edge).^2,2));
0266
0267 function O = omega(mdl)
0268 O = sqrt(3)*mdl.max_edge_length ./ (2*sqrt(2)*mdl.rad_circumsphere);
0269
0270 function NSR = tet_normalized_shape_ratio(mdl)
0271 NSR = 3 * (mdl.rad_insphere ./ mdl.rad_circumsphere);
0272
0273 function NSR = tri_normalized_shape_ratio(mdl)
0274 NSR = 2*mdl.rad_incircle./mdl.rad_circumcircle;
0275
0276 function A = dihedral_angles(mdl)
0277 v = nchoosek(1:4,2);
0278 for i = 1:6
0279 N1 = mdl.normals(mdl.elem2face(:,v(i,1)) , :) ...
0280 .* repmat(sign(mdl.inner_normal(:,v(i,1))-0.5),1,3);
0281 N2 = mdl.normals(mdl.elem2face(:,v(i,2)), :) ...
0282 .* repmat(sign(mdl.inner_normal(:,v(i,2))-0.5),1,3);
0283 C = cross3(N1, N2);
0284 D = -dot3 (N1, N2);
0285 A(:,i) = atan2( sqrt(sum(C.^2,2)), D);
0286 end
0287
0288
0289 function A = solid_angles(mdl)
0290 for i = 1:4
0291 E = mdl.elems';
0292 idx = 1:numel(E); idx(i:4:end) = [];
0293 N = mdl.nodes(E(idx),:) - reshape(repmat(mdl.nodes(E(i:4:end),:)',3,1),3,[])';
0294 nmrtr = abs(det3(N));
0295 L = sqrt(sum(N.^2,2));
0296 dnmtr = L(1:3:end).*L(2:3:end).*L(3:3:end) ...
0297 + dot3(N(1:3:end,:), N(2:3:end,:)) .* L(3:3:end) ...
0298 + dot3(N(1:3:end,:), N(3:3:end,:)) .* L(2:3:end) ...
0299 + dot3(N(2:3:end,:), N(3:3:end,:)) .* L(1:3:end);
0300
0301 A(:,i) = atan2( nmrtr, dnmtr );
0302 end
0303 idx = A<0;
0304 A(idx) = A(idx) + pi;
0305 A = 2*A;
0306
0307 function t = tri_angles(mdl)
0308 v = 1:3;
0309 L = mdl.edge_length(mdl.face2edge);
0310 L2 = L.^2;
0311 for i = 1:3
0312 t(:,i) = acos( (L2(:,v(1)) + L2(:,v(2)) - L2(:,v(3)) ) ...
0313 ./(2 .* L(:,v(1)) .* L(:,v(2)) ) );
0314 v = circshift(v,[1 2]);
0315 end
0316
0317
0318 function R = circumsphere_radius(mdl)
0319 ne = size(mdl.elems,1);
0320 E = mdl.elems';
0321 idx = 1:numel(E); idx(1:4:end) = [];
0322 N = mdl.nodes(E(idx),:) - reshape(repmat(mdl.nodes(E(1:4:end),:)',3,1),3,[])';
0323
0324
0325 X = repmat((1:3*ne)',1,3);
0326 T = reshape((1:3*ne)',3,[])';
0327 Y(1:3:3*ne,1:3) = T;
0328 Y(2:3:3*ne,1:3) = T;
0329 Y(3:3:3*ne,1:3) = T;
0330 A = sparse(X,Y,N);
0331 B = sum(N.^2,2);
0332 v = A\B;
0333 V = reshape(v,3,[])';
0334 R = 0.5 * sqrt(sum(V.^2,2));
0335 C = mdl.nodes(E(1:4:end),:) + 0.5*V;
0336
0337
0338 function H = tet_altitudes(mdl)
0339 elem_sorted = sort(mdl.elems,2);
0340 v = [4 3 2 1];
0341 for i = 1:4
0342 Nf = mdl.normals(mdl.elem2face(:,i),:);
0343 Cf = mdl.face_centre(mdl.elem2face(:,i),:);
0344 Pe = mdl.nodes(elem_sorted(:,v(i)), :);
0345 H(:,i) = abs(dot3(Nf, Pe - Cf));
0346 end
0347
0348 function H = tri_altitudes(mdl)
0349 s = sum(mdl.edge_length(mdl.face2edge),2) / 2;
0350 S = repmat(s,1,3);
0351 nmrtr = 2*sqrt(s .* prod(S - mdl.edge_length(mdl.face2edge),2));
0352 H = repmat(nmrtr,1,3) ./ mdl.edge_length(mdl.face2edge) ;
0353
0354
0355
0356 function R = circumcircle_radius(mdl)
0357 R = 0.25*prod(mdl.edge_length(mdl.face2edge),2) ./ mdl.face_area;
0358
0359 function R = incircle_radius(mdl)
0360 R = 2*mdl.face_area ./ sum(mdl.edge_length(mdl.face2edge),2);
0361
0362 function R = insphere_radius(mdl)
0363 R = 3* mdl.elem_vol ./ mdl.elem_area;
0364
0365 function A = elem_area(mdl)
0366 A = sum(mdl.face_area(mdl.elem2face),2);
0367
0368 function V = elem_volume(mdl)
0369 E = mdl.elems';
0370 idx = 1:numel(E); idx(1:4:end) = [];
0371 N = mdl.nodes(E(idx),:) - reshape(repmat(mdl.nodes(E(1:4:end),:)',3,1),3,[])';
0372 V = abs(det3(N))/6;
0373
0374 function do_unit_test
0375
0376 nodes = [0 0 0; 0 1 0; 1 1 0; 1 0 0;...
0377 0 0 1; 0 1 1; 1 1 1; 1 0 1];
0378 elems = [1 2 3 6; 3 6 7 8; 1 5 6 8; 1 3 4 8; 1 3 6 8];
0379 cube = eidors_obj('fwd_model','cube','nodes', nodes, 'elems', elems);
0380
0381 Q = calc_mesh_quality(cube);
0382
0383 if 0
0384 f = fieldnames(Q);
0385 for i = 1:length(f)
0386 disp(f{i});
0387 f2 = fieldnames(Q.(f{i}));
0388 for j = 1:length(f2)
0389 disp([' .' f2{j}]);
0390 disp(Q.(f{i}).(f2{j})');
0391 end
0392 end
0393 end
0394
0395 unit_test_cmp('CUBE:tri.NSR', Q.tri.NSR, 0.828427124746190, 1e-8);
0396 unit_test_cmp('CUBE:tri.mu', Q.tri.mu, 0.717438935214301, 1e-8);
0397 unit_test_cmp('CUBE:tri.eta', Q.tri.eta, 0.866025403784439, 1e-8);
0398 unit_test_cmp('CUBE:tri.theta', Q.tri.theta, 0.707106781186547, 1e-8)
0399 unit_test_cmp('CUBE:tri.iota', Q.tri.iota, 0.891518811420827, 1e-8)
0400 unit_test_cmp('CUBE:tri.kappa', Q.tri.kappa, 0.577350269189626, 1e-8)
0401 unit_test_cmp('CUBE:tri.min_angle', Q.tri.min_angle, 0.750000000000000, 1e-8)
0402 unit_test_cmp('CUBE:tet.NSR', Q.tet.NSR(1:4), 0.732050807568877, 1e-8);
0403 unit_test_cmp('CUBE:tet.mu', Q.tet.mu(1:4), 0.732050807568877, 1e-8);
0404 unit_test_cmp('CUBE:tet.eta', Q.tet.eta(1:4), 0.839947366596582, 1e-8);
0405 unit_test_cmp('CUBE:tet.min_angle', Q.tet.min_angle(1:4), 0.776074828029885, 1e-8);
0406
0407
0408 real = mk_common_model('b3t2r',[16,1]); real = real.fwd_model;
0409 [Q real] = calc_mesh_quality(real, 1);
0410 unit_test_cmp('B3T2:tri.NSR',Q.tri.NSR(1:3), [0.828427124746190;
0411 0.828051875000156; 0.828051875000156], 1e-8);
0412 unit_test_cmp('B3T2:tri.eta',Q.tri.eta(1:3), [0.866025403784439;
0413 0.865565849069265; 0.865565849069265], 1e-8);
0414 unit_test_cmp('B3T2:tet.min_angle',Q.tet.min_angle(1:3), [0.638037414014943;
0415 0.624804247692675; 0.638037414014943], 1e-8);
0416
0417 unit_test_cmp('B3T2:real.min_edge_length', real.min_edge_length(1:3), ...
0418 [19.190248174528644;18.575016823680134;18.575016823680134], 1e-8);
0419 unit_test_cmp('B3T2:real.rad_incircle', real.rad_incircle(7:9), ...
0420 [ 5.871326906131572; 5.440496467001706; 6.523833186814214], 1e-8);
0421
0422
0423 shell = real; shell.elems = shell.boundary;
0424 [Q shel] = calc_mesh_quality(shell, 1);
0425 unit_test_cmp('SHELL:tri.kappa', Q.tri.kappa(1:3), ...
0426 [ 0.577350269189626; 0.577043899379510; 0.577350269189626], 1e-8);
0427 unit_test_cmp('SHELL:dihedral_angle', shell.dihedral_angle(1:2,1:3), ...
0428 [ 0.830397390875299, 1.570796326794897, 1.570796326794897;
0429 1.570796326794897, 0.830397390875299, 1.570796326794897], 1e-8);
0430
0431 real = mk_library_model('pig_23kg_16el');
0432 [Q real] = calc_mesh_quality(real, 1);