function [Area] = Area_Q4_or_Q8 (n_n) % Revised 3/19/20 % Area of a Q4 or Q8 element by numerical integration % using natural coordinates interpolation and quadrature % number of nodes per element n_n = 4 or 8 only if (nargin == 0) ; n_n = 8 ; end ; % if assign a default n_e = 1 ;% number of elements in mesh (for plot) n_m = n_n ;% number of nodes in mesh (for plot) n_t = 1 ;% number of element types in mesh (for plot) % select which element model to use if ( n_n == 4 ) ; % number of nodes per element nodes = [1 2 3 4] ; % element connections x = [0. 48. 0. -60.] ; % x-coordinates y = [0. 64. 100. 40.] ; % y-coordinates n_q = 1 % number of quadrature pts % addpath /clear/www/htdocs/mech517/Akin_FEA_Lib % plot_input_2d_mesh (n_e, n_m, n_n, n_t, x, y, nodes); elseif ( n_n == 8 ) ; % alternate choice of quad nodes = [1 2 3 4 5 6 7 8];% connections x = [0. 48. 0. -60. 31. 29. -30. -30.] ;% x-coord y = [0. 64. 100. 40. 30. 86. 77. 31.] ;% y-coord n_q = 4 % number of pts % q8_mesh_plot (x, y, nodes) ; % plot curved el else ; % invalid user input error ('Only valid function input is 4 or 8') end ; % if straight or curved element picked % Set element coordinates as rectangular array Coord = [x; y]' ; % physical coordinates of nodes % get the quadrature data over (-1,-1) to (1,1) OK = [1, 4, 9, 16, 25] ; % the only valid n_q inputs if ( ~any (n_q == OK) ) ; % input not on required list error ('n_q invalid in Area_Q4_or_Q8.m') ; % stop end ; % if bad data switch n_q ; % branch on argument value case 1 ; % precision: exact for polynomial of degree 1 a_q = [0] ; b_q = [0] ; w_q = [4] ; case 4 ; % precision: exact for polynomial of degree 3 a_q = [-5.7735026918962573e-1 5.7735026918962573e-1 ... -5.7735026918962573e-1 5.7735026918962573e-1] ; b_q = [-5.7735026918962573e-1 -5.7735026918962573e-1 ... 5.7735026918962573e-1 5.7735026918962573e-1] ; w_q = [ 1.0000000000000000e+0 1.0000000000000000e+0 ... 1.0000000000000000e+0 1.0000000000000000e+0] ; case 9 ; % precision: exact for polynomial of degree 5 a_q = [-7.7459666924148340e-1 0.0000000000000000e+0 ... 7.7459666924148340e-1 -7.7459666924148340e-1 ... 0.0000000000000000e+0 7.7459666924148340e-1 ... -7.7459666924148340e-1 0.0000000000000000e+0 ... 7.7459666924148340e-1] ; b_q = [-7.7459666924148340e-1 -7.7459666924148340e-1 ... -7.7459666924148340e-1 0.0000000000000000e+0 ... 0.0000000000000000e+0 0.0000000000000000e+0 ... 7.7459666924148340e-1 7.7459666924148340e-1 ... 7.7459666924148340e-1] ; w_q = [3.0864197530864201e-1 4.9382716049382713e-1 ... 3.0864197530864201e-1 4.9382716049382713e-1 ... 7.9012345679012341e-1 4.9382716049382713e-1 ... 3.0864197530864201e-1 4.9382716049382713e-1 ... 3.0864197530864201e-1]; otherwise ; % use function that generated above list [a_q, b_q, w_q] = qp_rule_nat_quad (n_q); % use 1D data end ; % switch on number of points Area = 0.0 ; % initialiaze area sum for q = 1: n_q ; % loop over integration points ---> ---> a = a_q(q) ; b = b_q(q) ; w = w_q(q) ; % get pts & wt % Interpolate over (-1,-1) to (1,1) % H = element interpolation functions at a,b % DLH = interpolation local derivatives at a,b switch n_n % select from available library of elements case {4} % four node quadrilateral, -1 <= a,b <= 1, Q4 % (See qp_rule_nat_quad for corresponding data) % 4 3 % type interpolation functions, 1 x 4, 1-2-3-4 CCW % 1 2 H = [(1-a)*(1-b), (1+a)*(1-b), ... (1+a)*(1+b), (1-a)*(1+b)]/4 ; % element type parametric derivatives, 2 x 4 DLH = [-(1-b), (1-b), (1+b), -(1+b); ... % dH/da -(1-a), -(1+a), (1+a), (1-a)]/4 ; % dH/db case {8} % eight node quadrilateral, 1 x 8, b % (See qp_rule_nat_quad for tabulated data) % 4 7 3 % corners 1-2-3-4 CCW, mid 5-6-7-8. % 8 6 a % for-1 <= a,b <= 1 % 1 5 2 a_p = 1 + a; a_m = 1 - a; b_p = 1 + b; b_m = 1 - b; H = [a_m*b_m*(a_m+b_m-3), a_p*b_m*(a_p+b_m-3), ... a_p*b_p*(a_p+b_p-3), a_m*b_p*(a_m+b_p-3), ... 2*b_m*(1-a*a), 2*a_p*(1-b*b), ... 2*b_p*(1-a*a), 2*a_m*(1-b*b)]/4 ; % H(a,b) % element type parametric derivatives, 2 x 8 DLH = [-b_m*(a_m+a_m+b_m-3), b_m*(a_p+a_p+b_m-3), ... b_p*(a_p+a_p+b_p-3), -b_p*(a_m+a_m+b_p-3), ... -4*a*b_m, 2*(1-b*b), -4*a*b_p, -2*(1-b*b); % dH/da -a_m*(b_m+a_m+b_m-3), -a_p*(b_m+a_p+b_m-3), ... a_p*(b_p+a_p+b_p-3), a_m*(b_p+a_m+b_p-3), ... -2*(1-a*a), -4*b*a_p, 2*(1-a*a), -4*b*a_m]/4 ; % dH/db end ; % switch on number of nodes % calculate the Jacobian, for this point only Jacobian = DLH * Coord ; % 2 by 2 Jacobian matrix at q Det_J = det (Jacobian) ; % scalar determinate at q Area = Area + Det_J * w ; end ; % loop over integration points <--- <--- <--- <--- fprintf('Area of quadrilateral is %9.3e \n', Area)% result % end Area_Q4_or_Q8 ======================================= % Area of quadrilateral is 5.851e+03 function [r_q, s_q, w_q] = qp_rule_nat_quad (n_q) % =============== % tables of quadrature point locations and weights for % quadrilaterals interpolated in natural coords, -1 <= r,s <= 1 % n_1 = number of 1-D points in each direction % n_q = 2-D number of quadrature points required, = n_1^2 % r_q = all of first parametric quadrature coordinates % s_q = all of second parametric quadrature coordinates % w_q = weight at all quadrature points in n_p r_q = zeros (n_q, 1); s_q = zeros (n_q, 1) ; w_q = zeros (n_q, 1) ; if ( n_q == 1 ) ; % exact for constant or linear polynomial r_q = [ 0 ] ; s_q = [ 0 ] ; w_q = [ 4 ] ; % centroid point data % tensor products of 1-D rule elseif ( n_q == 4 | n_q == 9 | n_q == 16 | n_q == 25 ); % tabulated n_1 = fix ( sqrt ( n_q ) ) ; % size of 1-D rule [r_1, w_1] = qp_rule_Gauss (n_1) ; % get 1-D rule data k = 0 ; % initialize point number for i = 1 : n_1 ; % loop over points in s-direction for j = 1 : n_1 ; % loop over points in r-direction k = k + 1 ; % point number in the quadrilateral w_q (k) = w_1 (i) * w_1 (j) ; % product of 1-D weights r_q (k) = r_1 (j) ; % r-coordinate in quadrilateral s_q (k) = r_1 (i) ; % s-coordinate in quadrilateral end ; % for j point in quadrilateral end ; % for i point in quadrilateral else % update the tables or elseif above error ('\nERROR quad rule not tabulated for these points') end % if number of quadrature points % end qp_rule_nat_quad % ==========================================