# Example: Scordelis-Lo roof
#
# Loop over different discretizations
#
# -------------------------------------------------------------------

  elts = {"s3", "s4"};    % List of admissible element types

# Reference solution

  wref = -92.17;

# Discretization data

  etype = input("Enter element type (s3 or s4): ", "s");
  ne    = input("Enter array with number of elements: ");

  if (! ismember(etype, elts))
     error("Element type %s not supported\n", etype);
  end

# Output file

  fid = fopen("roof_loop.res", "wt");

  na = length(ne);

  for n = 1 : na

     fprintf(fid, "Discretization: ");
     fprintf(fid, "%d %s elements in one direction\n\n",
                  ne(n), etype);

# Model

     mdl = model(ne(n), etype);

# Create and export component

     roof = mfs_new(fid, mdl);

# Analysis

     roof = mfs_stiff(roof);
     roof = mfs_mass(roof, "consistent");
     roof = mfs_statresp(roof);

# Displacement at midpoint of free edge

     disp    = mfs_getresp(roof, "statresp", "disp", "G");
     wrel(n) = disp(3) / wref;

     clear mdl; clear roof;

  end

# Print convergence history

  fprintf(fid, "\nConvergence Data:\n\n");
  fprintf(fid, "  No. of Elements         wrel     \n")
  fprintf(fid, "  ---------------------------------\n");

  for n = 1 : na
      fprintf(fid, "     %5d              %7.4f\n", ne(n), wrel(n));
  end

# Plot convergence history

  figure(1, "position", [100, 100, 800, 500]);
  plot(ne, wrel, "marker", "o");
  grid;
  xlabel('Number of elements');
  ylabel('w_{rel}');

# Output convergence history

  data = [ne', wrel'];
  dlmwrite([etype, ".csv"], data, "precision", 5);

  fclose(fid);
