calc_mesh_quality

PURPOSE ^

CALC_MESH_QUALITY Various measures of mesh quality.

SYNOPSIS ^

function [Q mdl] = calc_mesh_quality(mdl, show)

DESCRIPTION ^

CALC_MESH_QUALITY Various measures of mesh quality.
  [Q MDL] = CALC_MESH_QUALITY(MDL) calculates several measures of mesh 
  quality for the EIDORS fwd_model struct MDL and returns them in Q 
  (described below).
  [Q MDL] = CALC_MESH_QUALITY(MDL, SHOW) if SHOW is TRUE displays two 
  figures with histograms of each measure.
  
   Q.tet.
     NSR    - inradius to circumradius ratio (normalized shape ratio)
                 n*r / R
                where r = inradius; R = circumradius; n = dimensions [3,4]
     mu     - inradius to longest edge ratio (L)                     [5]
                 r / L
     tau    - shortest to longest edge ratio                         [5]
                 l/L
     reg    - "regularity" of tetrahedron                      [6,p.151]
                 4r/H
                where H is the longest altitude
     zeta   - measure relating volume and area                       [7]
                 V^4 / ( SUM( Ai^2 ) ^3 )
                where V = volume; Ai = area of face i
     eta    - measure realting volume and edge length                [7]
                 V^(2/3) / SUM( li^2 )
                where li = length of edge i
     alpha  - ratio of volume to mean edge length                    [6]
                 V / (1/6 SUM( Ai ))
                where Ai = area of face i
     gamma  - ratio of volume to RMS edge length                     [8]
                 V / SQRT( 1/6 SUM( Ai^2 ) )
     min_angle - minimum dihedral angle

  Q.tri.
     NSR   - inradius to circumradius ratio (normalized shape ratio)
                 n*r / R
                where r = inradius; R = circumradius; n = dimensions [3,4]
     mu    - inradius to longest edge ratio (L)                        [9]
                 r / L       = tri_mu(mdl);
     eta   - measure relating triangle area and edge length        [10,11]
                 A / SUM( li^2) 
     theta - measure relating shortest altitude and edge length        [2]
                 h / SQRT( SUM ( li^2 ) )
     iota  - measure relating triangle area and edge length           [12]
                 A / SUM( li ) ^2
     kappa - ratio of shortest altitude to longest edge             [2,13]
                 h / L
     min_angle - minimum angle

  All the above measures are non-dimensional and scaled to be bounded 
  between 0 (worst) and 1 (best), where a regular tetrahedron/triangle
  scores 1.

  Measures where implemented based on review articles of V. Phathasarahy
  (1994) [1] and D. Field (2000) [2]. The references above indicate the
  original sources, although those were not consulted.

  Additionally, CALC_MESH_QUALITY returns the original MDL with several
  useful fields added.

  REFERENCES:
  [1] Parthasarathy V (1994) "A comparison of tetrahedron quality measures" 
      Finite Elem Anal Des 15:255–261
  [2] Field D (2000) "Qualitative measures for initial meshes" Int J Numer
      Meth Eng, 906:887–906
  [3] Cavendish JC, Field DA, Frey WH (1985) "An approach to automatic
      three-dimensional finite element mesh generation" Int J Numer Met 
      Eng 21:329-47
  [4] Field DA (1991) "A generic Delaunay triangulation algorithm for
      finite elemen meshes" Adv Eng Softw 13:263-72
  [5] Baker TJ (1989) "Element quality in tetrahedral meshes" Proc 7th Int
      Conf on Finite Element Methods in Flow Problems, Huntsville, AL.
      1018-24
  [6] Dannelongue HH, Tanguy PA (1991) "Three-dimensional adaptive finite 
      element computations and applications to non-newtonian flows" Int J 
      Numer Meth Fl 13:145–165
  [7] Cougny HL, Shephard MS, Georges MK (1990) "Explicit node point
      smoothing within Octree" Report No. 10-1990, SCOREC, RPI, Troy, NY.
  [8] Parthasarathy V, Kodiyalam (1991) "A constrained optimization
      approach to finite element mesh smoothing" Finite Elem Anal Des
      9:309-20
  [9] Perronnet A (1992) "Triangulation par arbre-4 de triangles equilateraux 
      at maximisation de la qualite" Publication du Laboratoire d’Analyse
      Numerique, Universite Pierr it Marie Curie et Centre National de la
      Recherche Scientifique
 [10] Bhatia RP, Lawrence KL (1990) "Two-dimensional finite element mesh 
      generation based on stripwise automatic triangulation" Comput Struct
      36:309–319.
 [11] Bank RE, Xu J (1996) "An algorithm for coarsening unstructured 
      meshes" Numer Math 73:1–36
 [12] Watabayshi GY, Galt JA (1986) "An optimized triangular mesh system 
      from random points" in "Numerical Grid Generation in Computational 
      Fluid Dynamics",Hauser J, Taylor C. (eds). Pineridge Press: Swansea,
      1986; 437–438
 [13] Suhara J, Fukuda J (1972) "Automatic mesh generation for finite 
      element analysis" on "Advances in Computational Methods in 
      Structural Mechanics and Design" Oden JT, Clough RW, Yamamoto Y (eds).
      University of Alabamal Press: Tuscaloosa, AL, 1972; 607–624

 See also FIX_MODEL

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [Q mdl] = calc_mesh_quality(mdl, show)
0002 %CALC_MESH_QUALITY Various measures of mesh quality.
0003 %  [Q MDL] = CALC_MESH_QUALITY(MDL) calculates several measures of mesh
0004 %  quality for the EIDORS fwd_model struct MDL and returns them in Q
0005 %  (described below).
0006 %  [Q MDL] = CALC_MESH_QUALITY(MDL, SHOW) if SHOW is TRUE displays two
0007 %  figures with histograms of each measure.
0008 %
0009 %   Q.tet.
0010 %     NSR    - inradius to circumradius ratio (normalized shape ratio)
0011 %                 n*r / R
0012 %                where r = inradius; R = circumradius; n = dimensions [3,4]
0013 %     mu     - inradius to longest edge ratio (L)                     [5]
0014 %                 r / L
0015 %     tau    - shortest to longest edge ratio                         [5]
0016 %                 l/L
0017 %     reg    - "regularity" of tetrahedron                      [6,p.151]
0018 %                 4r/H
0019 %                where H is the longest altitude
0020 %     zeta   - measure relating volume and area                       [7]
0021 %                 V^4 / ( SUM( Ai^2 ) ^3 )
0022 %                where V = volume; Ai = area of face i
0023 %     eta    - measure realting volume and edge length                [7]
0024 %                 V^(2/3) / SUM( li^2 )
0025 %                where li = length of edge i
0026 %     alpha  - ratio of volume to mean edge length                    [6]
0027 %                 V / (1/6 SUM( Ai ))
0028 %                where Ai = area of face i
0029 %     gamma  - ratio of volume to RMS edge length                     [8]
0030 %                 V / SQRT( 1/6 SUM( Ai^2 ) )
0031 %     min_angle - minimum dihedral angle
0032 %
0033 %  Q.tri.
0034 %     NSR   - inradius to circumradius ratio (normalized shape ratio)
0035 %                 n*r / R
0036 %                where r = inradius; R = circumradius; n = dimensions [3,4]
0037 %     mu    - inradius to longest edge ratio (L)                        [9]
0038 %                 r / L       = tri_mu(mdl);
0039 %     eta   - measure relating triangle area and edge length        [10,11]
0040 %                 A / SUM( li^2)
0041 %     theta - measure relating shortest altitude and edge length        [2]
0042 %                 h / SQRT( SUM ( li^2 ) )
0043 %     iota  - measure relating triangle area and edge length           [12]
0044 %                 A / SUM( li ) ^2
0045 %     kappa - ratio of shortest altitude to longest edge             [2,13]
0046 %                 h / L
0047 %     min_angle - minimum angle
0048 %
0049 %  All the above measures are non-dimensional and scaled to be bounded
0050 %  between 0 (worst) and 1 (best), where a regular tetrahedron/triangle
0051 %  scores 1.
0052 %
0053 %  Measures where implemented based on review articles of V. Phathasarahy
0054 %  (1994) [1] and D. Field (2000) [2]. The references above indicate the
0055 %  original sources, although those were not consulted.
0056 %
0057 %  Additionally, CALC_MESH_QUALITY returns the original MDL with several
0058 %  useful fields added.
0059 %
0060 %  REFERENCES:
0061 %  [1] Parthasarathy V (1994) "A comparison of tetrahedron quality measures"
0062 %      Finite Elem Anal Des 15:255–261
0063 %  [2] Field D (2000) "Qualitative measures for initial meshes" Int J Numer
0064 %      Meth Eng, 906:887–906
0065 %  [3] Cavendish JC, Field DA, Frey WH (1985) "An approach to automatic
0066 %      three-dimensional finite element mesh generation" Int J Numer Met
0067 %      Eng 21:329-47
0068 %  [4] Field DA (1991) "A generic Delaunay triangulation algorithm for
0069 %      finite elemen meshes" Adv Eng Softw 13:263-72
0070 %  [5] Baker TJ (1989) "Element quality in tetrahedral meshes" Proc 7th Int
0071 %      Conf on Finite Element Methods in Flow Problems, Huntsville, AL.
0072 %      1018-24
0073 %  [6] Dannelongue HH, Tanguy PA (1991) "Three-dimensional adaptive finite
0074 %      element computations and applications to non-newtonian flows" Int J
0075 %      Numer Meth Fl 13:145–165
0076 %  [7] Cougny HL, Shephard MS, Georges MK (1990) "Explicit node point
0077 %      smoothing within Octree" Report No. 10-1990, SCOREC, RPI, Troy, NY.
0078 %  [8] Parthasarathy V, Kodiyalam (1991) "A constrained optimization
0079 %      approach to finite element mesh smoothing" Finite Elem Anal Des
0080 %      9:309-20
0081 %  [9] Perronnet A (1992) "Triangulation par arbre-4 de triangles equilateraux
0082 %      at maximisation de la qualite" Publication du Laboratoire d’Analyse
0083 %      Numerique, Universite Pierr it Marie Curie et Centre National de la
0084 %      Recherche Scientifique
0085 % [10] Bhatia RP, Lawrence KL (1990) "Two-dimensional finite element mesh
0086 %      generation based on stripwise automatic triangulation" Comput Struct
0087 %      36:309–319.
0088 % [11] Bank RE, Xu J (1996) "An algorithm for coarsening unstructured
0089 %      meshes" Numer Math 73:1–36
0090 % [12] Watabayshi GY, Galt JA (1986) "An optimized triangular mesh system
0091 %      from random points" in "Numerical Grid Generation in Computational
0092 %      Fluid Dynamics",Hauser J, Taylor C. (eds). Pineridge Press: Swansea,
0093 %      1986; 437–438
0094 % [13] Suhara J, Fukuda J (1972) "Automatic mesh generation for finite
0095 %      element analysis" on "Advances in Computational Methods in
0096 %      Structural Mechanics and Design" Oden JT, Clough RW, Yamamoto Y (eds).
0097 %      University of Alabamal Press: Tuscaloosa, AL, 1972; 607–624
0098 %
0099 % See also FIX_MODEL
0100 
0101 % (C) 2013 Bartlomiej Grychtol. License: GPL v2 or v3.
0102 % $Id: calc_mesh_quality.m 5467 2017-05-10 10:39:02Z aadler $
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 % tet properties
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    % Note: The measure based on the minimum solid angle (sin(1/2 min angle)
0145    % is bound between 0 and 1, but both 0 and 1 are degenerate, so measure is
0146    % not calculated
0147    mdl.solid_angle    = solid_angles(mdl);
0148    mdl.dihedral_angle = dihedral_angles(mdl);
0149 end
0150 
0151 % tri properties
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 % TRI Quality measures
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 % only report measures for boundary triangles
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    % TET Quality measures
0175    Q.tet.NSR          = tet_normalized_shape_ratio(mdl);
0176    % this measure is not very sensitive and not bound between 0 and 1
0177    % Q.omega        = omega(mdl);
0178    Q.tet.mu           = mu(mdl);
0179    Q.tet.tau          = tau(mdl);
0180    Q.tet.reg          = R_measure(mdl);
0181    % These two measures from Field's paper are not scale invariant
0182    % Q.R_star       = R_star(mdl);
0183    % Q.S            = S_measure(mdl);
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 % subplot(3,3,9)
0212 % hist(180*real.dihedral_angle(:)/pi,100)
0213 % xlabel('dihedral angle');
0214 % xlim([0 180])
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); % choose 2 faces
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)); % length of each vector
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 % prepare a matrix for solving a system
0324 % A = sparse(length(N),length(N));
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; % center
0336 
0337 
0338 function H = tet_altitudes(mdl)
0339 elem_sorted = sort(mdl.elems,2);
0340 v = [4 3 2 1]; % which vertex is not on the face
0341 for i = 1:4 % loop over vertices
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 % note the the 5th element is a regular tet
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);
Generated on Tue 31-Dec-2019 17:03:26 by m2html © 2005