function []=Field_2D_Types_Ave_Lib(DEFAULT)% Revised 3/06/19 % DEFAULT omitted, require user interactive inputs % DEFAULT = 1 use quadratic triangle (T6) plane conductivity % Requires: addpath /clear/www/htdocs/mech517/Akin_FEA_Lib % -------------------------------------------------------- % 2D solution of the general field equation PDF % -[d/dx (Kx du/dx) + d/dy (Ky du/dy)] + G u - Q = 0 % for type 1 = conduction triangle T3 (or T6 or T10) % 2 = conduction quadrilateral Q4 (or Q8 or Q12) % 3 = flux NBC line element L2 (or L3 or L4) % 4 = mixed NBC line element L2 (or L3 or L4) % NOTE: If G=0 then this is the Poisson Equation % If G=0 and Q=0 then this is the Laplace Equation % If G=(u - u_inf), include G * u_inf in Q % -------------------------------------------------------- % TYPE 1 or 2 compatable area conduction elements: % With properties 1=Kx, 2=Ky, 3=G, 4=Q in el_prop(1-4) % and constitutive law: q_vector = +- [k] * Grad u % Iff Fourier = 1: use Darcy & Fourier laws q = -[k] * Grad u % else: Torsion etc q = +[k] * Grad u % Iff Q(x,y) source is defined at element nodes then % Q_2 = el_prop (5) ; % source at second node of the element % Q_3 = el_prop (6) ; % source at third node of the element % etc, so n_vals = n_n + 3 where % n_vals = number of property values for each element or type % if n_vals > 5 then Q is given at each element node % TYPE 3 edge elements with NBC known normal flux: % -Kn du/dn = q_in (postive when entering) % With properties 1=q_in, 2==0, 3==0, 4==0, to n_vals % TYPE 4 edge elements with NBC known mixed condition: % -Kn du/dn = g_in + h_in u % With properties 1=g_in, 2==h_in, 3==0, 4==0, to n_vals % NOTE: for convection g_in = h*u_infinity, h_in = h (conv coeff) %............................................................. % Mech 417, Finite Element Analysis, Rice University % Copyright J.E. Akin 2019. All rights reserved. %............................................................. % NOTE: input text files begin with msh_ and end with .txt % see file FEA_notation.m for variable descriptions % see file FEA_functions.m for functions summary % echo user information about this specific application echo_user_remarks % ===================================================================== % PHASE 1 % *** SET OR READ USER PROBLEM DEPENDENT CONTROLS *** % ===================================================================== if ( nargin == 0 ) % default to user defined field controls DEFAULT = 0 ; % require user inputs end % if use default application % set application controls if ( DEFAULT == 0 ) % then read user control inputs % axisym = 1 if the problem is axisymmetric % load_pt = 1 if point soures are input % n_g = number of DOF per node % N_N = number of EXPECTED nodes per element % n_p = dimension of parametric space % n_q = number of quadrature points required % n_r = number of rows in element B_e matrix % N_S = number of EXPECTED physical space dimensions [axisym, load_pt, n_g, N_N, n_p, n_q, n_r, N_S] = get_user_controls ; elseif ( DEFAULT == 1 ) % 2D T6 plane conductivity % set user maximum array sizes n_g = 1 ; % number of DOF per node N_N = 6 ; % number of EXPECTED nodes per element n_p = 2 ; % dimension of parametric space n_q = 4 ; % number of quadrature points n_r = 2 ; % number of rows in element B_e matrix N_S = 2 ; % EXPECTED physical space dimensions axisym = 0 ; % turn on(1)/off(0) the axisymmetric flag load_pt = 0 ; % turn on(1)/off(0) point source (load can be zero) end % if DEFAULT % set user logic choices (or add read function) color = 1 ; % turn on(1)/off(0) color plot of results as png file debug = 0 ; % turn on(1)/off(0) debug output listings echo_p = 1 ; % turn on(1)/off(0) listing of properties el_react = 0 ; % turn on(1)/off(0) listing of element reactions integral = 0 ; % turn on(1)/off(0) listing of solution integral Fourier = 1 ; % turn on(1)/off(0) minus sign in material law mass_ave = 0 ; % turn on(1)/off(0) ave of consistent and diagonal mass mass_dia = 0 ; % turn on(1)/off(0) use of scaled diagonal mass matrix mass_pt = 0 ; % turn on(1)/off(0) point mass (can be zero) post = 1 ; % turn on(1)/off(0) post-processing option show_bc = 0 ; % turn on(1)/off(0) list assembled S and C, after EBC show_e = 2 ; % turn on(>0)/off(0) list some S_e, M_e and C_e show_qp = 0 ; % turn on(1)/off(0) list qp post-processing results show_s = 0 ; % turn on(1)/off(0) listing of assembled S and C stiff_pt = 0 ; % turn on(1)/off(0) point stiffness (value can be zero) vector = 0 ; % turn on(1)/off(0) flux vector plot as png file % Is post-processing required ? if ( el_react > 0 | integral > 0 | show_qp > 0 | vector > 0 ) post = 1 ; % then post-processing is required end ; % if if ( Fourier == 0 ) % use plus sign law = 1 ; % torsion else % use minus sign law = -1 ; % Fourier Law or Darcy Law end ; % if sign of material law %% Echo all default or input controls %fprintf ('User control variables \n') %fprintf ('n_g number of DOF per node = %i \n', n_g) %fprintf ('n_n number of nodes per element = %i \n', N_N) %fprintf ('n_p dimension of parametric space = %i \n', n_p) %fprintf ('n_q number of quadrature points = %i \n', n_q) %fprintf ('n_r number of rows in B_e = %i \n', n_r) %fprintf ('n_s dimension of physical space = %i \n', N_S) %fprintf ('\n') % ===================================================================== % PHASE 2 % *** BEGIN DATA INPUT AND CHECKING PHASE *** % ===================================================================== fprintf ('Begin Phase 2, input mesh, properties, sources and EBC \n') % Read mesh nodal input data files % (from msh_bc_xyz.txt, with one line per node) % n_m = number of nodal points in mesh % n_s = number of space dimensions % PD = nodal packed digit essential BC flag (n_g digits 0 or 1) % NOTE: PD is either 0 or 1 here % x = mesh x-coordinate values % y = mesh y-coordinate values, if present, etc. for z [n_m, n_s, PD, x, y, z] = get_mesh_nodes (n_g) ; % bc_flags, coordinates % Compare actual data against space control if ( n_s ~= N_S ) % then possible control error fprintf ('WARNING: Read %i coordinates, but expected %i \n', n_s, N_S) end ; % if % *** Read mesh element types and nodal connections (padded with 0's) *** % (from msh_typ_nodes.txt, with one line per element) % n_e = number of elements in the mesh % n_n = maximum number of nodes per element % el_type = type number of every element, usually 1 % nodes = node connection list for each element (pad with 0 to size n_n) % n_t = number of different element types [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements ; % Compare actual data against connection control if ( n_n ~= N_N ) % number of expected nodes fprintf ('WARNING Read %i nodes per element but expected %i \n', n_n, N_N) end ; % if inconsistent inputs % Allocate element type control data if more than one type input type_i = zeros (n_t, 1) ; % maximum DOF per el type type_n = zeros (n_t, 1) ; % maximum nodes per el type type_p = zeros (n_t, 1) ; % dimension of type parametric space type_q = zeros (n_t, 1) ; % number of type quadrature points type_x = zeros (n_t, 1) ; % number of type geometry nodes (default type_n) % Default to a single element type if ( n_t == 1 ) % then initialize type_n = n_n ; type_p = n_p ; type_q = n_q ; type_x = n_n ; else % read controls for each type ! NEW ******** NEW % *** Read each element type control integers *** % (from msh_ctrl_types.txt, with one line per element type) % type_n = maximum nodes per el type % type_p = dimension of type parametric space % type_q = number of type quadrature points % type_x = number of type geometry nodes (default type_n) [type_n, type_p, type_q, type_x] = get_el_type_controls (n_t) ; % check type control vs system control has_n = max ( type_n ) ; % actual max node per el has_p = max ( type_p ) ; % actual max parm dim has_q = max ( type_q ) ; % actual max Gauss rule has_x = max ( type_x ) ; % actual max geom nodes if ( has_n > n_n | has_x > n_n ) % warn fprintf ('Warning: type_n = %i is inconsistent with mesh controls \n', has_n) n_n = has_n % increase the maximum number of nodes per element end ; % if if ( has_p > n_p ) % warn fprintf ('Warning: type_p = %i is inconsistent with mesh controls \n', has_p) n_p = has_p % increase the maximum parametric space dimension end ; % if if ( has_q > n_q ) % warn fprintf ('Warning: type_q = %i is inconsistent with mesh controls \n', has_q) n_q = has_q % increase the maximum number of quadrature points end ; % if end ; % if more than one element type type_i = n_g * type_n ; % number of degrees of freedom per type % Read any material properties for the elements % (from msh_properties.txt, one line per element or type) % n_mats = number of materials, = n_e, or n_t if homogeneous % n_vals = number of property values for each element or type % msh_properties (n_mats, n_vals) = system properties array [n_mats, n_vals, msh_properties] = get_mesh_properties (echo_p, n_e, n_t) ; fprintf ('Above properties are: \n') fprintf ('TYPE 1 or 2: Kx, Ky, G, Q (Q2,...,Qn) \n') fprintf ('For -[d/dx (Kx du/dx) + d/dy (Ky du/dy)] + G u - Q = 0 \n') fprintf ('TYPE 3: 1=q_in, 2==0, 3==0, 4==0 \n') fprintf ('For -Kn du/dn = q_in (postive when entering) \n') fprintf ('TYPE 4: 1=g_in, 2==h_in, 3==0, 4==0 \n') fprintf ('For -Kn du/dn = g_in + h_in * u, where for convection: \n') fprintf ('h_in = h the convection coeff), g_in = h * u_infinity \n') % Count maximum local and system equation numbers n_d = n_g * n_m ; % system degrees of freedom (DOF) n_i = n_g * n_n ; % number of DOF per element % Pre-allocate and initialize system matrices and vectors. % Ans = solution vector % C = total source vector % C_pt = point source part of C for reaction recovery % EBC_flag = n_g individual binary flags at each node % EBC_value = EBC value for each restrained DOF % el_prop = element or type properties % M = mass matric (if any) % S = stiffness matrix % Type_integral = integral of solution (G*u) in each element % Type_i_G_u = integral of scaled solution (G*u) in each element Ans = zeros (n_d, 1) ; % results C = zeros (n_d, 1) ; % force or source C_pt = zeros (n_d, 1) ; % point source portion of C C_react = zeros (n_d, 1) ; % copy of C for reaction recovery EBC_flag = zeros (n_m, n_g) ; % initialize all DOF flags EBC_value = zeros (n_m, n_g) ; % default EBCs to zero el_prop = zeros (n_vals, 1); % initialize element or type properties M = zeros (n_d, n_d) ; % generalized mass Q_v = zeros (n_n, 1) ; % variable source nodal values S = zeros (n_d, n_d) ; % generalized stiffness Type_integral = zeros (n_t, 1) ; % integral of the solution Type_i_G_u = zeros (n_t, 1) ; % integral of G times the solution Type_volume = zeros (n_t, 1) ; % volume of each type of element % Read point sources & couples data, if any, and insert in C % (from file msh_load_pt.txt, one line per point source) if ( load_pt > 0 ) ; % need point sources data [C_pt] = get_point_sources (n_g, n_m) ; % get point sources if ( axisym == 1 ) % then source acts over radius for i = 1:n_m rows = [1:n_g] + n_g*(i-1) ; % get rows at point C_pt (rows) = C_pt (rows) * x (i) ; % times radius end ; % for end ; % if C = C + C_pt ; % add point sources end ; % if any point source expected % Read point stiffness data, if any, and add to S diagonal % (from file msh_stiff_pt.txt, one line per point stiffness) if ( stiff_pt > 0 ) ; % need point stiffness data [S] = get_and_add_point_stiffness (n_g, n_m, S); % get point stiffness end ; % if any point stiffness expected % Extract and count EBC flags from packed integer flag PD [EBC_flag] = get_ebc_flags (n_g, n_m, PD) ; % unpack flags % EBC_count = number of essential (Dirichlet) BC's EBC_count = sum( sum ( EBC_flag == 1 ) ) ; % number of EBC if ( EBC_count == 0 ) fprintf ('WARNING: No essential boundary condition given. \n') else fprintf ('Note: expecting %i essential BC values. \n', EBC_count) end ; % if % Read essential boundary condition (EBC) values, if any % (from file msh_ebc.txt one line per constraint) if ( EBC_count > 0 ) ; % need EBC data [EBC_value] = get_ebc_values (n_g, n_m, EBC_flag) ; % read data else % unusual, but possible fprintf ('WARNING: no essential boundary conditions flagged \n') end ; % if any EBC data expected % Prepare a default general 2D mesh plot as png file plot_input_2d_mesh (n_e, n_m, n_n, n_t, x, y, nodes) pause (2) % wait 2 seconds before next plot % Set an integer to simplify element type selection below domain = n_n * 100 + n_g * 10 + n_p ; % parametric reference code % *** END OF DATA INPUT AND CHECKING PHASE *** % echo current controls and logic as debugging aid if ( debug == 1 ) % provide debug checks final_control_status (n_d, n_e, n_g, n_i, n_m, n_mats, n_n, ... n_p, n_q, n_r, n_s, n_t, n_vals) final_logic_status (axisym, color, echo_p, el_react, integral, ... Fourier, load_pt, mass_ave, mass_dia, mass_pt, post, ... show_bc, show_e, show_qp, show_s, stiff_pt, vector) else fprintf ('NOTE: no debug outputs requested. \n') end ; % debug listings % ====================================================================== % PHASE 3 % ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM ** % ====================================================================== fprintf ('Begin Phase 3, assemble (scatter) element arrays \n') if ( show_e > 0 ) % then list some element matrices fprintf ('Displaying the first %i sets of element matrices \n', show_e) end ; % if arrays to display % Open unit needed for post-processing operations if ( post == 1 ) % open post-processing element unit el_unit = fopen ('el_store.bin', 'w') ; % open for element writing if ( el_unit == -1 ) % open failed fprintf ('\nWARNING: post-processing unit open failed \n') post = 0 ; % turn off post-processing option end ; % if unit available qp_unit = fopen ('qp_store.bin', 'w') ; % open for GQ writing if ( qp_unit == -1 ) % open failed fprintf ('\nWARNING: post-processing unit open failed \n') post = 0 % turn off post-processing option end ; % if unit available end ; % if open post-processing unit type_old = 0 ; % initialize for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> type = el_type (j) ; % current element type number typ_i = type_i (type) ; % maximum DOF per el type typ_n = type_n (type) ; % maximum nodes per el type typ_p = type_p (type) ; % dimension of type parametric space typ_q = type_q (type) ; % number of type quadrature points domain = typ_n * 100 + n_g * 10 + typ_p ; % parametric reference code e_nodes = nodes (j, 1:typ_n) ; % connectivity % Gather nodal coordinates xy_e (1:typ_n, 1) = x(e_nodes(1:typ_n)) ; % x coord at el nodes if ( n_s > 1 ) % then 2D or 3D xy_e (1:typ_n, 2) = y(e_nodes(1:typ_n)) ; % y coord at el nodes end ; % if 2D or 1D % Allocate and clear element type arrays C_e = zeros (typ_i, 1) ; % clear array el sources C_r = zeros (typ_i, 1) ; % clear array el reactions E_e = zeros (n_r, n_r) ; % material constitutive matrix H_i = zeros (typ_n, 1) ; % scalar interpolations integral Jacobian = zeros (n_s, n_s) ; % geometry Jacobian (one-to-one) M_e = zeros (typ_i, typ_i) ; % clear array el convection M_g = zeros (typ_i, typ_i) ; % clear generalized mass matrix S_e = zeros (typ_i, typ_i) ; % clear array el stiffness % Set constant element type properties if ( n_mats > n_t ) % each element has some different properties J = j ; % the element number else % use type data J = type ; % the type number end ; % if which property values el_prop (1:n_vals) = msh_properties (J, 1:n_vals) ; % Shorthand type property names switch ( type ) % branch on element type case {1, 2} % compatable continuum elements for PDE % -[d/dx (Kx du/dx) + d/dy (Ky du/dy)] + G u - Q = 0 % u = unknown % Kx = el_prop (1) ; % modulus % Ky = el_prop (2) ; % modulus % G = el_prop (3) ; % convection or foundation like term % Q = el_prop (4) ; % internal source term % n_vals = number of property values for each element or type % if n_vals > 5 then addition sources are: % Q_2 = el_prop (5) ; % source at second node of the element % Q_3 = el_prop (6) ; % source at third node of the element, etc if ( n_vals < 4 ) % then property data error error ('This application requires 4 properties') end ; % if n_vals Kx = el_prop (1) ; % modulus Ky = el_prop (2) ; % modulus if ( Kx == 0 | Ky == 0 ) % then invalid data error ('Kx and/or Ky is zero. Invalid properties.') end ; % if G = el_prop (3) ; % convection or foundation like term Q = el_prop (4) ; % internal source term Q_v (1:n_n) = Q ; % if constant source if ( n_vals == 3+n_n ) % then variable Q(x, y) at each elem node for Q_k = 1:n_n Q_v (Q_k) = el_prop (Q_k+3) ; % source at each node end ; % for end ; % if loading % Fill constitutive array E_e = [ Kx 0 ; 0 Ky] ; % orthotropic case {3} % edge element with known source % TYPE 3 edge elements with knowmn flux: % -Kn du/dn = q_in (+ when entering) % With properties 1=q_in, 2==0, 3==0, 4==0 q_in = el_prop (1) ; % nornal flux case {4} % edge element with mixed (Robin) condition % TYPE 4 edge elements with mixed condition: % -Kn du/dn = g_in + h_in u % With properties 1=g_in, 2==h_in, 3==0, 4==0 g_in = el_prop (1) ; % usually convection coeff * ref temp h_in = el_prop (2) ; % usually convection coeff otherwise error ('Type number of element exceeds contol value \n') end ; % switch on type number % Allocate and clear quadrature point arrays B_q = zeros (n_r, typ_i) ; % initalize B matrix values, at q DGH_q = zeros (n_s, typ_n) ; % interpolation global derivatives, at q DLH_q = zeros (typ_p, typ_n) ; % interpolation local derivatives, at q Jacobian = zeros (typ_p, n_s) ; % geometry Jacobian H_q = zeros (1, typ_n) ; % clear scalar interpolations, at q xy_q = zeros (1, n_s) ; % physical location of point q % ELEMENT, FOUNDATION STIFFNESSES AND SOURCE RESULTANT MATRICES % Get quadrature rule for current element type if ( type ~= type_old ) % then [r_q, s_q, t_q, w_q] = get_quadrature_rule (typ_n, typ_p, typ_q) ; end % if % Numerical Intergation of S_e, C_e , in 1D parametric volume = 0 ; % element volume for q = 1:typ_q ; % begin quadrature loop ---> ---> ---> ---> ---> ---> r = r_q (q) ; w = w_q (q) ; % recover data s = s_q (q) ; t = 0 ; % not used in 1D % Element scalar interpolation, H, & local derivatives, DLH [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t) ; % Interpolate global position of quadrature point xy_q = H_q * xy_e (1:typ_n, :) ; % 1 x n_s % global (x, y) % Compute the geometry mapping data, d_parm to d_physical [J_det, Jacobian] = parametric_to_physical (n_s, typ_n, typ_p, DLH_q, xy_e); if ( axisym ) % then axisymmetric-unit-radian surface J_det = xy_q (1) * J_det ; % theorem of Papus end ; % if axisymmetric body if ( type == 1 | type == 2 ) % then continuum element J_Inv = inv (Jacobian) ; % n_s x n_s % inverse % Global derivatives of H are DGH DGH_q = J_Inv * DLH_q ; % n_s x n_n % dH/dx & dH/dy % Insert gradient terms into the differential operator [B_q] = DGH_q ; % same for scalars % Update element square matrices S_e = S_e + (B_q' * E_e * B_q) * J_det * w_q (q) ; % conduction % Generalized mass matrix update (for variable Q) M_g = M_g + (H_q' * H_q) * J_det * w_q (q) ; % convection % Constant source C_e = C_e + H_q' * Q * J_det * w_q (q) ; % source % Update integral of the interpolation functions, for post-process H_i = H_i + H_q' * J_det * w_q (q) ; % H integral volume = volume + J_det * w_q (q) ; % element volume % save data needed for post-processing stresses if ( post == 1 ) % save data for this point fwrite (qp_unit, xy_q, 'double') ; % save coordinates fwrite (qp_unit, E_e, 'double') ; % save material fwrite (qp_unit, B_q, 'double') ; % save strain-displacement end ; % if save point data elseif ( type == 3 ) % then a known normal flux % -K_n * U,n + g_in = 0 C_e = C_e + H_q' * q_in * J_det * w_q (q) ; elseif ( type == 4 ) % then a mixed (Robin) condition % -K_n * U,n + h_in * U + g_in = 0 S_e = S_e + (H_q' * h_in * H_q) * J_det * w_q (q) ; C_e = C_e + H_q' * g_in * J_det * w_q (q) ; else error ('Element arrays not defined for type > 3 \n') end ; % if element type end ; % for quadratures <--- <--- <--- <--- <--- <--- <--- <--- <--- Type_volume (type) = Type_volume (type) + volume ; if ( type < 3 & n_vals == 3+n_n ) % then overwrite for variable source C_e = M_g * Q_v (1:n_n) ; % revised source resultant end ; % if if ( show_e > 0 ) % then list some element matrices if ( j <= show_e ) % show the first few fprintf ('Element matrix S_e = \n') ; disp(S_e) fprintf ('Element matrix C_e = \n') ; disp(C_e) end ; % if elements wanted end ; % if show_e if ( G ~= 0 ) % then convection or elastic foundation M_e = G * M_g ; % true surrounding matrix S_e = S_e + M_e ; % net sq matrix if ( show_e > 0 ) % then list some element matrices if ( j <= show_e ) % show the first few fprintf ('Element matrix M_e = \n') ; disp(M_e) fprintf ('Element S_e + M_e = \n') ; disp(S_e) end ; % if elements wanted end ; % if show_e end ; % if mass % Save element arrays needed for post-processing stresses if ( post == 1 ) % save arrays for this element fwrite (el_unit, xy_e, 'double') ; % save coordinates fwrite (el_unit, el_prop, 'double') ; % save properties fwrite (el_unit, S_e, 'double') ; % save beam stiffness fwrite (el_unit, M_g, 'double') ; % save generalized mass fwrite (el_unit, M_e, 'double') ; % save foundation stiffness fwrite (el_unit, C_e, 'double') ; % save line load resultant fwrite (el_unit, H_i, 'double') ; % save interpolation integral end ; % if save element arrays % SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS % Insert completed element matrices into system matrices [rows] = get_element_index (n_g, typ_n, e_nodes); % eq numbers % rows = vector subscript converting element to system eq numbers S (rows, rows) = S (rows, rows) + S_e ; % add to system stiffness C (rows) = C (rows) + C_e ; % add to sys force/couples % *** END SCATTER TO SYSTEM ARRAYS *** end ; % for each j element in mesh <<==== <<==== <<==== <<==== <<==== if ( show_s == 1 ) % then list each system matrix fprintf ('System matrix S = \n') ; disp(S) fprintf ('System matrix C = \n') ; disp(C) end ; % if show_s fprintf ('Total input source terms = %g \n', sum (C)) % *** END ELEMENT GENERATION % ======================================================================== % PHASE 4 % APPLY BC AND CONSTRAINTS (Make problem unique) % ======================================================================== fprintf ('Begin Phase 4, modify matrices for EBC \n') % ALLOCATE STORAGE FOR OPTIONAL REACTION RECOVERY if ( EBC_count > 0 ) ; % reactions occur % EBC_row = row from S for each EBC, for later reaction calculation % EBC_col = row from C for each EBC, for later reaction calculation EBC_row = zeros (EBC_count, n_d) ; % reaction data EBC_col = zeros (EBC_count, 1) ; % reaction data [EBC_row, EBC_col] = save_reaction_matrices (EBC_flag, S, C); end ; % if essential BC exist (almost always true) % in the future, inforce all MPC equations here % ENFORCE ESSENTIAL BOUNDARY CONDITIONS (and insert identities) if ( EBC_count > 0 ) ; % reactions occur [S, C] = enforce_essential_BC (EBC_flag, EBC_value, S, C); end ; % if essential BC exist (almost always true) if ( show_bc == 1 ) % then list each system matrix after EBC fprintf ('After EBC system matrix S = \n') ; disp(S) fprintf ('After EBC system matrix C = \n') ; disp(C) end ; % if show_s % ======================================================================== % PHASE 5 % SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS, OR % SOLVE EIGENVALUE SYSTEM % ======================================================================== fprintf ('Begin Phase 5, solve modified system \n') Ans = S \ C ; % Compute nodal unknowns list_save_results (n_g, n_m, Ans) ; % save for plots and print %color = 0 ; % skip link to plot library if ( color == 1 ) % display and save a png plot %addpath /clear/www/htdocs/mech517/Matlab_Plots % library color_scalar_result (x, y, nodes, Ans) ; % plot pause (2) ; % pause for 2 seconds scalar_result_surface (x, y, nodes, Ans) ; % plot pause (2) ; % pause for 2 seconds end ; % if plot desired % OPTIONAL FORCE REACTION RECOVERY & SAVE if ( EBC_count > 0 ) ; % reactions exist ? % EBC_react = nodal source required to maintain EBC value EBC_react = zeros (EBC_count, 1) ; % pre-allocate [EBC_react] = recover_reactions_print_save (n_g, n_d, ... EBC_flag, EBC_row, EBC_col, Ans); % reaction to EBC end ; % if EBC existed % ===================================================================== % PHASE 6 % POST-PROCESS ELEMENTS % Get location and values of derivaties at all q-point, and/or % Get integral of the solution and/or % Sum all the element or nodal values % ===================================================================== fprintf ('Begin Phase 6, use the answers \n') if ( post == 1 ) % then post-process all element stresses if ( qp_unit > 0 ) ; % close GQ written files fclose (qp_unit) ; % file qp_store.bin qp_unit = fopen ('qp_store.bin', 'r') ; % open for reading end ; % if if ( el_unit > 0 ) ; % close written binary files fclose (el_unit) ; % file el_store.bin el_unit = fopen ('el_store.bin', 'r') ; % open for reading end ; % if element arrays were saved if ( el_unit == -1 ) % open failed fprintf ('\nWARNING: post-processing unit open failed \n') end ; % if unit available if ( qp_unit == -1 ) % open failed fprintf ('\nWARNING: post-processing unit open failed \n') end ; % if unit available % Save temp file for flux plot input (see Matlab quiver_* plot files) q_id = fopen ('el_qp_xyz_fluxes.txt', 'w') ; % flux plot file fprintf ('\n') fprintf ('Element Post-processing: \n') strain = zeros (n_r, 1) ; % pre-allocate strains stress = zeros (n_r, 1) ; % pre-allocate stresses % Re-compute or read arrays from storage format short e for j = 1:n_e ; % loop over elements ====>> ====>> ====>> ====>> ====>> type = el_type (j) ; % current element type number typ_i = type_i (type) ; % maximum DOF per el type typ_n = type_n (type) ; % maximum nodes per el type typ_q = type_q (type) ; % number of type quadrature points if ( type > 2 ) % then skip this element and go to the next one fprintf ('NOTE: skipped type %i element number %i \n', type, j) continue % to the next j ====>> ====>> ====>> end ; % skip this type % Gather answers e_nodes = nodes (j, 1:typ_n) ; % connectivity [rows] = get_element_index (n_g, typ_n, e_nodes) ; % eq numbers Ans_e (1:typ_i) = Ans(rows) ; % gather element results % Recover previously calculated element arrays [xy_e] = fread (el_unit, [typ_n, n_s], 'double'); % coordinates [el_prop] = fread (el_unit, [n_vals, 1], 'double'); % properties [S_e] = fread (el_unit, [typ_i, typ_i], 'double'); % beam stiffness [M_g] = fread (el_unit, [typ_i, typ_i], 'double'); % generalized mass [M_e] = fread (el_unit, [typ_i, typ_i], 'double'); % foundation stiffness [C_e] = fread (el_unit, [typ_i, 1], 'double'); % line load resultant [H_i] = fread (el_unit, [typ_n, 1], 'double'); % interpolation integral for q = 1:typ_q ; % begin quadrature loop ---> ---> ---> ---> ---> ---> % Recover previous calculations at this point (or re-compute) [xy_q] =fread(qp_unit, [1, n_s], 'double'); % get coordinates [E_e] =fread(qp_unit, [n_r, n_r], 'double'); % get material matrix [B_q] =fread(qp_unit, [n_r, typ_i], 'double'); % get application def % Compute the generalized strains & stresses strain = B_q * Ans_e(1:typ_i)' ; % in-plane strains (gradient) stress = law * E_e * strain ; % in-plane stresses (fluxes) if ( axisym ) % then axisymetric solid stress = stress / xy_q (1) ; end ; % if % List quadrature point results fprintf ('\n') fprintf ('El, Q_Pt, Coordinates %i, %i, %g %g \n', ... j, q, xy_q (1:2)) fprintf ('El, Q_Pt, Gradient Vector %i, %i, %g %g \n', ... j, q, strain (1:2)) fprintf ('El, Q_Pt, Flux Vector %i, %i, %g %g \n', ... j, q, stress(1:2)) % save flux for plot programs fprintf (q_id, '%g %g %g %g \n', xy_q (1:2), stress (1:2)); end ; % for q quadrature points <--- <--- <--- <--- <--- <--- <--- <--- if ( integral == 1 & n_g == 1 ) % report solution integral Type_integral (type) = Type_integral (type) + dot (Ans_e, H_i) ; Type_i_G_u (type) = Type_i_G_u (type) + dot (Ans_e, H_i)*G ; end ; % if scalar solution % Optional output of nodal reactions on every element if ( el_react ) % list nodal heat flows for this element % Assign any point source to the first element with that node if ( load_pt ) % then point sources to assign (once) for eq = 1:n_i % loop over element nodes k = rows (eq) ; % find system DOF number if ( C_pt (k) ~= 0 ) % found an input point value C_e (eq) = C_e (eq) + C_pt (k); % copy to this node C_pt (k) = 0 ; % prevent a second use end ; % if first appearance of this point source in the system end ; % for eq node of the current element end ; % if load_pt has been assigned %b fprintf ('\nNode List for Element %i: \n', j) ; disp (e_nodes) %b fprintf ('Input Resultant Nodal Sources \n') ; disp (C_e') % Finally, get the element reactions and list them C_m = S_e * Ans_e' - C_e ; fprintf ('Reactions at Nodes \n') ; disp (C_m') end ; % if el_react list end ; % for each j element in mesh <<==== <<==== <<==== <<==== <<==== % close the flux data plot file fclose(q_id) %vector = 0 ; % skip link to plot library if ( vector == 1 ) % plot flux vectors at selected quadrature points %addpath /clear/www/htdocs/mech517/Matlab_Plots % library quiver_qp_flux_mesh (1.0, 1, law) % inputs (scale, inc_g, sign) pause (4) ; % pause for 4 seconds end ; % if if ( axisym ) % then axisymmetric solid Type_integral = 2 * pi() * Type_integral ; Type_i_G_u = 2 * pi() * Type_i_G_u ; Type_volume = 2 * pi() * Type_volume ; end ; % if if ( integral == 1 & n_g == 1) % report solution integral fprintf ('\n') if ( n_t > 1 ) % report on each type for j = 1:n_t fprintf ('Solution integral in type %i elements = %g \n', ... j, Type_integral (j)) fprintf ('Scaled integral in type %i elements = %g \n', ... j, Type_i_G_u (j)) fprintf ('Solution volume in type %i elements = %g \n', ... j, Type_volume (j)) end ; % for each type end ; % if multiple element types fprintf ('Total solution integral = %g \n', sum (Type_integral)) fprintf ('Total scaled integral = %g \n', sum (Type_i_G_u )) end ; % if scalar solution fprintf ('\n Total domain volume = %g \n', sum (Type_volume)) format short % OPTIONAL POST PROCESSING COMPLETED else fprintf ('\nNo element post-processing requested \n') end ; % for post-processing stresses % ===================================================================== % PHASE 7 Average nodal flux values % ( Also needed as input for an error estimate) % ===================================================================== %B free up memory used by system S & C matrices fprintf ('Begin SCP averaging flux at the nodes \n') % build database of element neighbors needed for node averages % n_e can be reduced above, use latest [L_first, L_last] = first_last_L_at_pt (n_e, nodes) ; % appear [L_to_L_sum] = count_elems_at_elem (n_e, L_first, L_last, nodes); % count [L_to_L_set] = form_elems_at_el (n_e, L_first, L_last, ... % set of nodes, L_to_L_sum) ; % neighbors % use patches of neighbors to average flux components at mesh nodes [scp_ave] = calc_ave_mesh_fluxes (n_e, n_q, n_r, n_s, nodes, x, y, z, ... L_to_L_sum, L_to_L_set) ; % get extreme range of each flux component [range, n_range] = scp_ave_flux_max_min (scp_ave) ; % amount & where ave_test = 2 * eps ; % small value if ( max (max (abs (range)) < ave_test )) ; % then warn fprintf ('WARNING: all averages are zero. \n') fprintf (' Error estimator will be omitted. \n') end % if warn % list average flux values if desired, save for error estimates show_ave = 1 ; % control name if ( show_ave == 1 ) ; % list ave fluxes list_ave_fluxes (scp_ave) ; % at each node end ; % if % save average flux values for plots and/or error estimates save_ave_fluxes (scp_ave) ; % to scp_node_ave_fluxes.txt % plot flux vectors in file scp_node_ave_fluxes.txt quiver_scp_flux_mesh ; % plot averaged flux vectors % ===================================================================== % PHASE 8 Estimate errors from flux averages % ===================================================================== % end % function Field_2D_Types_Ave_Lib % =============================