function [Area] = Area_T6_curved (n_q);% Revised 3/17/20 % Area of a curved T6 element by numerical integration % using unit coordinates interpolation and quadrature % n_q = number of quadrature points required % Area = integral dx dy (= integral 1 dx dy) % Area = integral |J(r,s)| dr ds % If straight with uniform nodes then |J| degree = 0 % If curved or non-uniform spaces |J| degree = 1+1 = 2 % % Area = 0 + Sum_(q=1, n_q) |J(r_q, s_q)| * w_q % For degree = 0, 1, 2, 3, 4, 5, 6, ... % need n_q >= 1, 1, 3, 4, 6, 7, 12, ... n_n = 6 ;% number of nodes per element 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) % Define vertices then mid-sides, counter-clockwise (why) nodes = [1 2 3 4 5 6] ; % connections x = [0. 48. 0. 31. 29. -20.] ; % x-coord y = [0. 64. 100. 30. 86. 50.] ; % y-coord % Set element coordinates as rectangular array Coord = [x; y]' ; % physical coordinates of nodes % plot the curved triangle as a data check % addpath /clear/www/htdocs/mech517/Akin_FEA_Lib % t6_mesh_plot (x, y, nodes) ; % plot curved T6 % get the quadrature data over (0,0) to (1,0) to (0,1) OK = [1, 3, 4] ; % the only valid n_q inputs here if ( ~any (n_q == OK) ) ; % input not on required list fprintf ('Non-existant rule %i requested \n', n_q) error ('n_q invalid in function Area_T6_curved') end ; % if bad data % Allocate storage for quadrature data, and extract them 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 = [ 1 / 3.] ; s_q = [ 1 / 3.] ; w_q = [ 1 / 2.] ; elseif ( n_q == 3 ) % exact for 2-rd degree polynomials r_q=[6.6666666666666667e-01, 1.6666666666666667e-01, ... 1.6666666666666667e-01 ] ; s_q=[1.6666666666666667e-01, 6.6666666666666667e-01, ... 1.6666666666666667e-01 ] ; w_q=[1.6666666666666666e-01, 1.6666666666666667e-01, ... 1.6666666666666667e-01 ] ; elseif ( n_q == 4 ) % exact for 3-rd degree polynomials r_q = [ 1.0, 0.6, 0.6, 1.8 ] / 3. ; s_q = [ 1.0, 0.6, 1.8, 0.6 ] / 3. ; w_q = [ -27, 25, 25, 25 ] / 96. ; end ; % if which rule to extract Area = 0.0 ; % initialiaze area sum for q = 1: n_q ; % loop over integration points ---> ---> r = r_q(q) ; s = s_q(q) ; w = w_q(q) ; % get pts & wt % H = element interpolation functions at r, s % DLH = interpolation local derivatives at r, s H = [(1-3*r-3*s+2*r*r+4*r*s+2*s*s), (2*r*r - r), ... (2*s*s-s), 4*(r-r*r-r*s), 4*r*s, 4*(s-r*s-s*s)]; DLH = [(-3+4*r+4*s), (-1+4*r), 0, ... (4-8*r-4*s), (4*s), (-4*s) ; ... % dH/dr (-3+4*r+4*s), 0, (-1+4*s), ... (-4*r), (4*r), (4-4*r-8*s)] ; % dH/ds % 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 ; % increment area value end ; % loop over integration points <--- <--- <--- <--- fprintf('Quadrilateral area is %10.4e \n', Area);% result % end Area_T6_curved ====================================