# Example: Scordelis-Lo roof
#
# Detailed analysis of one discretization
#
# --------------------------------------------------------------

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

# Reference solution

  wref = -92.17;

# Output device depending on operating system

  if (isunix)
     DEVICE = "svg";
  else
     DEVICE = "png";
  endif
  FORMAT = ["-d", DEVICE];
  EXT    = [".", DEVICE];

# Discretization data

  etype = input("Enter element type (s3 or s4): ", "s");
  nelt  = input("Enter number of elements in one direction: ");

  if (! ismember(etype, elts))
     error("Element type %s not supported\n", etype);
  end
  if (nelt <= 0)
     error("Number of elements must be strictly positive\n");
  end

# Output file

  fid = fopen([etype, ".res"], "wt");

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

# Model

  mdl  = model(nelt, etype);

# Create and export component

  roof = mfs_new(fid, mdl);
  mfs_export([etype, ".msh"], "msh", roof, 
             "mesh", "mesh", "axes");

# Analysis

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

# Displacements of free edge FG

  mfs_print(fid, roof, "statresp", {"disp", "FG"});

  dsp  = mfs_getresp(roof, "statresp", "disp", "G");
  wrel = dsp(3) / wref;
  fprintf(fid, "\nNormalized displacement: uz = %7.4f\n", wrel);

# Export displacements

  mfs_export([etype, ".dsp"], "msh", roof, "statresp", "disp");

# Element results

  roof = mfs_results(roof, "statresp", "element");
  mfs_print(fid, roof, "statresp", {"stress", "FG"}, 
            {"resultant", "FG"});
  mfs_export([etype, ".sig"], "msh", roof, "statresp", "stress");

# Stress resultants along edge y = L

  rslt = mfs_getresp(roof, "statresp", "resultant", "FG");
  r    = cell2mat(rslt);
  coor = reshape([r.coor], 3, length(rslt));
  phi  = 180 * atan2(coor(1, :), coor(3, :)) / pi; 

  titletxt = sprintf("%s elements, %d intervals",
                     etype, nelt);

  figure(1, "position", [100, 100, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  plot(phi, [r.Nx], "color", "green",
       phi, [r.Ny], "color", "red");
  title(titletxt);
  legend("N_x", "N_y", "location", "northwest");
  grid;
  xlabel('\phi [°]'); ylabel('[N/mm]');

  print(["N", EXT], FORMAT);

  figure(2, "position", [200, 100, 750, 500],
            "paperposition", [0, 0, 15, 10]);
  plot(phi, [r.Mx] / 1000, "color", "green",
       phi, [r.My] / 1000, "color", "red");
  title(titletxt);
  legend("M_x", "M_y", "location", "northeast");
  grid;
  xlabel('\phi [°]'); ylabel('[Nm/mm]');

  print(["M", EXT], FORMAT);

  fclose(fid);
