function Parametric_Line   % Revised 1/18/2012  
% Extract physical geometric data for a parametric line mapped
%        to straight line segments on the x-axis 
%   or to straight line or curved segments the xy-plane

%............................................................. 
%     Mech 417, Finite Element Analysis, Rice University 
%     Copyright J.E. Akin 2012. All rights reserved.
%............................................................. 
% NOTE: input text files begin with msh_ and end with .tmp
% see file FEA_notation.m for variable descriptions
global axisym  ; % the axisymmetric mesh flag, global variable
 
% echo information about input and output text files
% echo_files_info   % uncomment column 1, if desired

% echo user information about this specific application
echo_user_remarks

% =====================================================================
%                             PHASE 1
%           *** SET OR READ USER PROBLEM DEPENDENT CONTROLS *** 
% =====================================================================
%   set default user numerical values
n_g = 1 ; % number of degrees of freedom per node (fake here)
N_N = 2 ; % number of EXPECTED nodes per element
n_p = 1 ; % dimension of parametric space
n_q = 5 ; %4 ; %3 ; % number of quadrature points required
%           Exact for polynomial of degree (2n_q - 1)
N_S = 1 ; % number of EXPECTED physical space dimensions

%   set user logic choices
axisym = 0  ; % turn on(1)/off(0) an axisymmetric mesh flag

%    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_s dimension of physical space   = %i \n', N_S)
if ( axisym == 1 )
  fprintf ('The physical space is axisymmetric. \n')
end % if

% =====================================================================
%                               PHASE 2
%              *** BEGIN DATA INPUT AND CHECKING PHASE ***
% =====================================================================
%    Read mesh nodal input data files 
%    (from msh_bc_xyz.tmp, 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 fake (un-used) in this application, so enter 0
% x    = mesh x-coordinate values
% y    = mesh y-coordinate values, if present, etc. for z
%   Compare actual data against controls
[n_m, n_s, PD, x, y, z] = get_mesh_nodes (n_g) ; % bc_flags, coordinates 
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 ***
%   (from msh_typ_nodes.tmp, 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
% n_t     = number of different element types
[n_e, n_n, n_t, el_type, nodes] = get_mesh_elements  ; 
%   Compare actual data against controls
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

%   Check integration rule for moment of inertia
if ( n_n == 2 & n_q < 1 ) % increase quadrature points
  fprintf ('NOTE increased n_q to 1 \n')
  n_q = 1 ;
end % if a linear element
if ( n_n == 3 & n_q < 3 ) % increase quadrature points
  fprintf ('NOTE increased n_q to 3 \n')
  n_q = 3 ;
end % if a quadratic element
if ( n_n == 4 & n_q < 5 ) % increase quadrature points
  fprintf ('NOTE increased n_q to 5 \n')
  n_q = 5 ;
end % if a cubic element

%   Prepare a default general 2D mesh plot (uncomment as desired)
%   using /clear/www/htdocs/mech517/Matlab_Plots/mesh_shrink_plot.m
% mesh_shrink_plot % show a space between elements
% pause (3) % wait 3 seconds before next plot
%   or using /clear/www/htdocs/mech517/Matlab_Plots/mesh_plot.m
% mesh_plot    % no space between elements, mesh_plot (inc_e, inc_p)
%   or using true element shape with
%      /clear/www/htdocs/mech517/Matlab_Plots/L2_mesh_plot.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/L3_mesh_plot.m
%   or /clear/www/htdocs/mech517/Matlab_Plots/L4_mesh_plot.m

%  Set an integer to simplify element type selection below
domain = n_p * 100 + n_g * 10 + n_n ; % parametric reference code

%    *** END OF DATA INPUT AND CHECKING PHASE ***

% ======================================================================
%                              PHASE 3
%       ** GENERATE ELEMENT VALUES AND ASSYMBLE INTO SYSTEM **
% ======================================================================

%    *** Initially fill the quadrature data files ***
% r_q, s_q, t_q = parametric coordinates of quadrature point
% w_q           = parametric weight      of quadrature point
if ( n_q > 0 )                 % numerical integration to be used
  r_q = zeros (n_q, 1) ; s_q = zeros (n_q, 1) ;    % allocate
  t_q = zeros (n_q, 1) ; w_q = zeros (n_q, 1) ;    % allocate
  if ( n_p == 1 )                   % line elements
    s_q = 0 ; t_q = 0 ;
%   comment out one of the next two lines, re interpolation coords
    [r_q, w_q] = qp_rule_unit_Gauss (n_q) ; % if using unit coords
%   [r_q, w_q] = qp_rule_Gauss (n_q)      ; % if natural    coords
  elseif ( n_p == 2 )               % triangles or quadrilaterals
    error ('\n ERROR fill 2-D quadrature rules')
  else
    error ('\n ERROR fill 3-D quadrature rules')
  end % if parametric space dimension
end % if numerical integration data required 
 
%     Allocate and clear element arrays
  H_q   = zeros (n_n, 1)      ; % clear element interpolations
  DLH_q = zeros (n_s, n_n)    ; % interpolation local derivatives 

%     Initilize global integrals
Total_length = 0 ; Total_first_x = 0 ; Total_sec_yy = 0 ; 
Total_first_y = 0 ; Total_sec_xx = 0 ; Total_sec_xy = 0 ;
 
%                *** BEGIN ELEMENT LOOPS ***
for j = 1:n_e  ; % loop over elements ====>> ====>> ====>> ====>> 
  e_nodes = nodes (j, 1:n_n)                      ; % connectivity 
  type    = el_type (j)            ; % current element type number
   
%       Gather nodal coordinates (using vector subscripts)
  xy_e (1:n_n, 1) = x(e_nodes)             ; % x coord at el nodes 
  xy_e (1:n_n, 2) = y(e_nodes)             ; % y coord at el nodes 

%     Numerical integration of along line
  length = 0 ; first_x = 0 ; sec_yy = 0 ; % initialize summation
  first_y = 0 ; sec_xx = 0 ; sec_xy = 0 ;

%           *** BEGIN QUADRATURE LOOP IN CURRENT ELEMENT ***
  for q = 1:n_q   ; % begin quadrature loop ---> ---> ---> ---> ---> --->
    r = r_q (q) ; s = 0 ; t = 0 ;  w = w_q (q)     ; % recover data

%     Evaluate element interpolation data at (r, s, t), via domain code
    [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t) ;

%     Global position and differential lengths (Jacobian), at point q
    xy_q  = H_q * xy_e                      ;  % global (x, y)
    dx_dr = DLH_q * xy_e (:, 1)             ;  % dx (r) / dr  
    dy_dr = DLH_q * xy_e (:, 2)             ;  % dy (r) / dr  

    if ( n_p == 1 ) % then a parametric line
      if ( n_s == 1 ) % then a one-to-one map
        J_det = dx_dr                       ;  % dL (r) / dr
      elseif ( n_s == 2 ) % then a planar curve
        J_det = sqrt (dx_dr^2 + dy_dr^2)    ;  % dL (r) / dr
      else % need equation for space curve Jacobian
        error ('Add 3D line Jacobian calculations')
      end % if n_s, physical space dimension
    else % parametric surface or volume, n_p = 2 or 3
      error ('Missing parametric surface or volume Jacobian')
    end % if n_p, dimension of parametric space

    if ( axisym ) % then axisymmetric-unit-radian surface
      J_det = xy_q (1) * J_det ; % theorem of Papus
    end % if axisymmetric body 

    if ( J_det <= 0 ) % then inverse mapping not possible
      fprint ('Fatal mesh distortion at element %i \n', j)
      error  ('Zero or negative geometric Jacobian %g \n', J_det)
    end % if fatally distorted element 

%    Properties at this point, if spatially variable (none to get here)

%    Update integrals of the interpolation functions, for post-process
    length  = length  + J_det * w_q (q)                 ; % length
    first_x = first_x + xy_q (1)   * J_det * w_q (q)    ; % first moment
    first_y = first_y + xy_q (2)   * J_det * w_q (q)    ; % first moment
    sec_yy  = sec_yy  + xy_q (1)^2 * J_det * w_q (q)    ; % second moment
    sec_xx  = sec_xx  + xy_q (2)^2 * J_det * w_q (q)    ; % second moment
    sec_xy  = sec_xy  - xy_q (1)*xy_q (2)*J_det*w_q (q) ; % second moment

   end % for quadratures   <--- <--- <--- <--- <--- <--- <--- <--- <--- 
%         *** END OF QUADRATURE LOOP IN CURRENT ELEMENT ***

   if ( length <= 0 ) % then bad data
     fprintf ('Check connectivity data for element %i \n', j)
     error   ('Element has a negative length')
   end % if length

%         SCATTER TO (ASSEMBLE INTO) SYSTEM ARRAYS 

    Total_length  = Total_length  + length   ; % update system sums
    Total_first_x = Total_first_x + first_x  ; % update system sums
    Total_first_y = Total_first_y + first_y  ; % update system sums
    Total_sec_yy  = Total_sec_yy  + sec_yy   ; % update system sums
    Total_sec_xx  = Total_sec_xx  + sec_xx   ; % update system sums
    Total_sec_xy  = Total_sec_xy  + sec_xy   ; % update system sums

%       *** END SCATTER TO SYSTEM ARRAYS ***
 
end % for each j element in mesh   <<==== <<==== <<==== <<==== <<==== 

%       *** END ELEMENT GENERATION, Post process

% ========================================================================
%                          PHASE 4
%         APPLY BC AND CONSTRAINTS (Make problem unique)
% ========================================================================
%      None for this application

% ========================================================================
%                          PHASE 5
%         SOLVE LINEAR EQUATION SYSTEM AND ITS REACTIONS, OR
%                  SOLVE EIGENVALUE SYSTEM
% ========================================================================
%      None for this application

% =====================================================================
%                             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
% =====================================================================
%      find the centroid from system sums
x_bar = Total_first_x / Total_length ;
y_bar = Total_first_y / Total_length ;

%      Output the system results
fprintf ('\n')
fprintf ('Straight line results: \n')
fprintf ('Length           = %g \n', Total_length)
fprintf ('X-centroid point = %g \n', x_bar)
fprintf ('Global I_yy      = %g \n', Total_sec_yy)
if (n_s > 1 ) % then more results
  fprintf ('Y-centroid point = %g \n', y_bar)
  fprintf ('Global I_xx      = %g \n', Total_sec_xx)
  fprintf ('Global I_xy      = %g \n', Total_sec_xy)
end % if n_s

% =====================================================================
%                             PHASE 7
%           ESTIMATE ERROR IN EACH ELEMENT AND THE SYSTEM
% =====================================================================
%      None for this application

% end function Parametric_Line   % =======================
   
%         ++++++++++++++++++++++++++++++++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
%  ++++++++ active functions in alphabetical order ++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
%         ++++++++++++++++++++++++++++++++++++++++++
 
% See /clear/www/htdocs/mech517/Matlab_Plots for graphic options 
% and http://www.clear.rice.edu/mech517/help_plot.html for help files
 
function echo_files_info  % ==========================
file_id = fopen('msh_files_info.tmp', 'r'); % open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_files_info.tmp) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf (remark_line)        ; % echo the line
  end % while data exists
end % if file supplied
fprintf('\n')
% end function echo_file_info  % ==========================

function echo_user_remarks  % ==========================
file_id = fopen('msh_remarks.tmp', 'r'); ;% open for read only
if ( file_id > 0 )                 % then file exists
  fprintf('\n(Echo of msh_remarks.tmp) \n')
  while feof(file_id) ==0          % read until an end of file
    remark_line = fgets(file_id) ; % read a string and append \n
    fprintf (remark_line)        ; % echo the line
  end % while data exists
end % if file supplied
fprintf('\n')
% end function echo_user_remarks  % ==========================

function [H_q, DLH_q] = el_shape_n_local_deriv (domain, r, s, t)
% parametric interpolation functions and their derivatives
% at local point r, or at r, s, or at r, s, t, and 
% domain = n_p * 100 + n_g * 10 + n_n ; % an integer code
% Currently using unit coordinates, 0 <= r,s,t <= 1 

 switch domain     % select from available library of elements
% One-dimensional elements
  case {112}       % two node line linear element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q = [(1 - r),  r] ; DLH_q = [-1,    1] ;

  case {113}       % three node line quadratic element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q   = [(1 - 3*r + 2*r^2), (4*r - 4*r^2), (2*r^2 - r)] ;
    DLH_q = [(-3 + 4 * r),      (4 - 8*r),     (4*r -1)] ;

  case {114}       % four node line cubic element
%   (See qp_rule_unit_Gauss for corresponding tabulated data)
    H_q   = [(1 - r*11/2 + 9*r^2 - 9*r^3/2), (9*r - 45*r^2/2 + 27*r^3/2), ...
             (-9*r/2 + 18*r^2 - 27*r^3/2),   (r - 9*r^2/2 + 9*r^3/2) ] ;
    DLH_q = [(-11/2 + 18*r - 27*r^2/2),      (9 - 45*r + 81*r^2/2), ...
             (-9/2 + 36*r - 81*r^2/2),       (1 - 9*r + 27*r^2/2) ] ;
  otherwise
    fprintf ('Current mesh domain flag = %i \n', domain)
    error ('\n Missing parametric functions for that domain')
 end % switch
% end function el_shape_n_local_deriv 

function [n_e, n_n, n_t, el_type, nodes] = get_mesh_elements (); 
% input file controls (for various data generators)  
 
load  msh_typ_nodes.tmp        ; % el_type, connectivity list (3)  
n_e = size (msh_typ_nodes,1)   ; % number of elements  
if ( n_e == 0 )                ; % data file missing  
  error ('\n Error missing file msh_typ_nodes.tmp')  
end % if error  
n_n = size (msh_typ_nodes,2) - 1 ; % nodes per element        
fprintf ('\n')
fprintf ('Read %i elements with type & %i nodes each. \n', ...
          n_e, n_n)  
fprintf ('(Echo element number, file msh_typ_nodes.tmp) \n')
% Extract all element type numbers into vector el_type
el_type = round (msh_typ_nodes(:, 1)); % el type number >= 1  
% Extract all element connection lists into array nodes
nodes (1:n_e, 1:n_n) = msh_typ_nodes (1:n_e, 2:1+n_n);  

if ( n_n > 9 ) % skip element number display
 disp(msh_typ_nodes (:, 1:1+n_n)) % echo data , without el number
else
 % Echo data with pretty print from control values
 for i = 1:n_e % loop over all elements ===> ===> ===> ===> ===>
  if ( n_n == 2 ) % nodes per element
   fprintf ('%i, %i %i %i       \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 3 )
   fprintf ('%i, %i %i %i %i    \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 4 )
   fprintf ('%i, %i %i %i %i %i \n', i, el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 6 )
   fprintf ('%i, %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 8 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  elseif ( n_n == 9 )
   fprintf ('%i, %i %i %i %i %i %i %i %i %i %i \n', ...
             el_type(i), nodes(i,1:n_n))
  end % if number of nodes to show
 end % for all i elements <=== <=== <=== <=== <=== <=== <=== <=== 
end % if show element number

n_t  = max(el_type)                   ; % number of element types  
if ( n_t > 1 )
  fprintf ('Note: maximum element type is %i \n', n_t)
end % if
% end get_mesh_elements    % ==========================
 
function [n_m, n_s, PD, x, y, z] = get_mesh_nodes (n_g); % =======
% input file controls (for various data generators)  
  
%    READ MESH AND EBC_FLAG SEQUENTIAL INPUT DATA  
load msh_bc_xyz.tmp       ; % bc_flag, x-, y-, z-coords  
n_m = size (msh_bc_xyz,1) ; % number of nodal points in mesh  
if ( n_m == 0 )           ; % data missing !  
  error ('\n Error missing nodal file msh_bc_xyz.tmp')  
end % if error  

n_s = size (msh_bc_xyz,2) - 1 ; % number of space dimensions  
fprintf ('Read %i nodes. \n', n_m)  
fprintf ('(Echo node, file msh_bc_xyz.tmp) \n')
fprintf ('node, bc_flags, %i coordinates \n', n_s)

msh_bc_xyz (:, 1)= round (msh_bc_xyz (:, 1)); % require integer 
PD = msh_bc_xyz (1:n_m, 1) ; % integer Packed BC flag  
x  = msh_bc_xyz (1:n_m, 2) ; % extract x column   
y (1:n_m, 1) = 0.0 ;  z (1:n_m, 1) = 0.0 ; % default to zero 
  
if (n_s > 1 )                         ; % check 2D or 3D 
  y = msh_bc_xyz (1:n_m, 3)           ; % extract y column   
end % if 2D or 3D 
if ( n_s == 3 )                       ; % check 3D 
  z = msh_bc_xyz (1:n_m, 4)           ; % extract z column  
end % if 3D  

% Echo data with pretty print from control values
for j = 1:n_m  % loop over all nodes ===> ===> ===> ===> ===> ===>
  if ( n_s == 1 )
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g \n', j, PD(j), x(j))
    else
      fprintf ('%i, %i    %g \n', j, PD(j), x(j))
    end % if dof per node
  elseif (n_s == 2)
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g %g \n', j, PD(j), x(j), y(j))
    else
      fprintf ('%i, %i    %g %g \n', j, PD(j), x(j), y(j))
    end % if dof per node
  else % n_s = 3
    if ( n_g == 2 )
      fprintf ('%i, %2.2i %g %g %g \n', j, PD(j), x(j), y(j), z(j))
    else
      fprintf ('%i, %i    %g %g %g \n', j, PD(j), x(j), y(j), z(j))
    end % if dof per node
  end % if space dimension
end % loop for j-th node <=== <=== <=== <=== <=== <=== <=== <===
% end get_mesh_nodes    % ==========================
 
function [r_q, w_q] = qp_rule_Gauss (n_q) % ================
%  tables of quadrature point locations and weights
%  for lines interpolated in natural coordinates
%             -1 <= r <= 1,   sum (w_q) == 2
 switch n_q
  case 1 % precision: exact for polynomial of degree 1
    r_q (1) = 0 ; w_q = 2 ;

  case 2 % precision: exact for polynomial of degree 3
    r_q (1) = -1 / sqrt (3) ; r_q (2) = 1 / sqrt (3) ;
    w_q (1) = 1             ; w_q (2) = 1 ;

  case 3 % precision: exact for polynomial of degree 5
    r_q (1) = -0.774596669241483377035835 ;
    r_q (2) =  0.000000000000000000000000 ;
    r_q (3) =  0.774596669241483377035835 ;
    w_q (1) =  0.55555555555555555555556 ;
    w_q (2) =  0.88888888888888888888889 ;
    w_q (3) =  0.55555555555555555555556 ;

  case 4 % precision: exact for polynomial of degree 7
    r_q (1) = -0.861136311594052575223946 ;
    r_q (2) = -0.339981043584856264802666 ;
    r_q (3) =  0.339981043584856264802666 ;
    r_q (4) =  0.861136311594052575223946 ;
    w_q (1) =  0.34785484513745385737306 ;
    w_q (2) =  0.65214515486254614262694 ;
    w_q (3) =  0.65214515486254614262694 ;
    w_q (4) =  0.34785484513745385737306 ;

  case 5 % precision: exact for polynomial of degree 9
    r_q (1) = -0.906179845938663992797627 ;
    r_q (2) = -0.538469310105683091036314 ;
    r_q (3) =  0.000000000000000000000000 ;
    r_q (4) =  0.538469310105683091036314 ;
    r_q (5) =  0.906179845938663992797627 ;
    w_q (1) =  0.23692688505618908751426 ;
    w_q (2) =  0.47862867049936646804129 ;
    w_q (3) =  0.56888888888888888888889 ;
    w_q (4) =  0.47862867049936646804129 ;
    w_q (5) =  0.23692688505618908751426 ;

  otherwise % update the tables
   fprintf ('See text for tabulated data where n_q = %i \n', n_q)
   error ('\n ERROR quadrature rule not tabulated for these points')
 end % switch number of quadrature points
% end qp_rule_Gauss % ===============================================

function [r_q, w_q] = qp_rule_unit_Gauss (n_q) % ================
%  tables of quadrature point locations and weights
%  for lines interpolated in unit coordinates. 
%              0 <= r <= 1,   sum (w_q) == 1
 switch n_q
  case 1 % precision: exact for polynomial of degree 1
    r_q (1) = 0.5 ; w_q = 1 ;

  case 2 % precision: exact for polynomial of degree 3
    r_q (1) = 2.1132486540518711774543e-01 ;
    r_q (2) = 7.8867513459481288225457e-01 ;
    w_q (1) = 5.0000000000000000000000e-01 ;
    w_q (2) = 5.0000000000000000000000e-01 ;

  case 3 % precision: exact for polynomial of degree 5
    r_q (1) = 1.1270166537925831148208e-01 ;
    r_q (2) = 5.0000000000000000000000e-01 ;
    r_q (3) = 8.8729833462074168851792e-01 ;
    w_q (1) = 2.7777777777777777777778e-01 ;
    w_q (2) = 4.4444444444444444444445e-01 ;
    w_q (3) = 2.7777777777777777777778e-01 ;

  case 4 % precision: exact for polynomial of degree 7
    r_q (1) = 6.9431844202973712388027e-02 ;
    r_q (2) = 3.3000947820757186759867e-01 ;
    r_q (3) = 6.6999052179242813240133e-01 ;
    r_q (4) = 9.3056815579702628761197e-01 ;
    w_q (1) = 1.7392742256872692868653e-01 ;
    w_q (2) = 3.2607257743127307131347e-01 ;
    w_q (3) = 3.2607257743127307131347e-01 ;
    w_q (4) = 1.7392742256872692868653e-01 ;

  case 5 % precision: exact for polynomial of degree 9
    r_q (1) = 4.6910077030668003601186e-02 ;
    r_q (2) = 2.3076534494715845448184e-01 ;
    r_q (3) = 5.0000000000000000000000e-01 ;
    r_q (4) = 7.6923465505284154551816e-01 ;
    r_q (5) = 9.5308992296933199639881e-01 ;
    w_q (1) = 1.1846344252809454375713e-01 ;
    w_q (2) = 2.3931433524968323402065e-01 ;
    w_q (3) = 2.8444444444444444444444e-01 ;
    w_q (4) = 2.3931433524968323402065e-01 ;
    w_q (5) = 1.1846344252809454375713e-01 ;

  otherwise % update the tables
   fprintf ('See text for tabulated data where n_q = %i \n', n_q)
   error ('\n ERROR quadrature rule not tabulated for these points')
 end % switch number of quadrature points
% end qp_rule_unit_Gauss % ===============================================

%                       END SOURCE FILES
%                          EXAMPLE

%  >> Parametric_1d_line

%  (Echo of msh_remarks.tmp) 
%  ============== Begin Application Remarks ====================
%  Numerical evaluation of the geometric properties of a line
%  from x=2 to x=6 so length is 4, cg is 4, inertia is 69.33
%  Nodal dof: none (geometry study only)
%  Element type = 1 (a parametric line)
%  Element connection: two nodes per element
%  Element properties: none (later add areas)
%  NOTE: nodal BC_Flag not used here, enter zero for it
%  =============== End Application Remarks =======================

%  User control variables 
%  n_g number of DOF per node        = 1 
%  n_n number of nodes per element   = 2 
%  n_p dimension of parametric space = 1 
%  n_q number of quadrature points   = 3 
%  n_s dimension of physical space   = 1 

%  Read 5 nodes. 
%  (Echo node, file msh_bc_xyz.tmp) 
%  node, bc_flags, 1 coordinates 
%  1, 0  2  0     
%  2, 0  3  0     
%  3, 0  4  0     
%  4, 0  5  0     
%  5, 0  6  0     

%  Read 4 elements with type & 2 nodes each. 
%  (Echo element number, file msh_typ_nodes.tmp) 
%  1, 1 1 2 
%  2, 1 2 3 
%  3, 1 3 4 
%  4, 1 4 5 

%  Straight line results: 
%  Length           = 4 
%  X-centroid point = 4 
%  Global I_xx      = 69.3333 

%==============================================================
%  User control variables 
%  n_g number of DOF per node        = 1 
%  n_n number of nodes per element   = 2 
%  n_p dimension of parametric space = 1 
%  n_q number of quadrature points   = 3 
%  n_s dimension of physical space   = 1 

%  Read 5 nodes. 
%  (Echo node, file msh_bc_xyz.tmp) 
%  node, bc_flags, 2 coordinates 
%  1, 0    2 1.5 
%  2, 0    3 2.25 
%  3, 0    4 3 
%  4, 0    5 3.75 
%  5, 0    6 4.5 

%  Read 4 elements with type & 2 nodes each. 
%  (Echo element number, file msh_typ_nodes.tmp) 
%  1, 1 1 2       
%  2, 1 1 2       
%  3, 1 1 2       
%  4, 1 1 2       

%  WARNING: Read 2 coordinates, but expected 1 

%  Read 5 mesh coordinate pairs 
%  Read 4 elements with 2 nodes each 

%  Straight line results: 
%  Length           = 5 
%  X-centroid point = 4 
%  Y-centroid point = 3 
%  Global I_yy      = 86.6667 
%  Global I_xx      = 48.75 
%  Global I_xy      = -65 
%==============================================================