# Example: Scordelis-Lo roof
#
# Detailed analysis of one discretization and comparison of results
# obtained with different element types
#
# -------------------------------------------------------------------

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

# Reference solution

  wref = -92.17;

# Discretization data

  ne = input("Enter number of elements in one direction: ");

  if (ne <= 0)
     error("Number of elements must be strictly positve\n");
  end

# Output file

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

# Loop over element types

  for n = 1 : length(elts)

      etype = elts{n};

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

      mdl  = model(ne, etype);

      roof = mfs_new(fid, mdl);

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

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

      roof = mfs_results(roof, "statresp", "element");

      switch etype

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

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

      endswitch

      clear roof; clear mdl;

  endfor

# Plot results

  figure(1, "position", [100, 100, 1000, 759],
            "paperposition", [0, 0, 15, 10]);
  plot(phi3, [s3.Nx], "color", "green", "marker", "^",
       phi3, [s3.Ny], "color", "red", "marker", "^",
       phi4, [s4.Nx], "color", "green", "marker", "s",
       phi4, [s4.Ny], "color", "red", "marker", "s");
  legend("N_x", "N_y", "N_x", "N_y", "location", "northwest");
  grid;
  xlabel('\phi [°]'); ylabel('[N/mm]');

  print("N.svg", "-dsvg");

  figure(2, "position", [200, 200, 1000, 759],
            "paperposition", [0, 0, 15, 10]);
  plot(phi3, [s3.Mx] * 1e-3, "color", "green", "marker", "^",
       phi3, [s3.My] * 1e-3, "color", "red", "marker", "^",
       phi4, [s4.Mx] * 1e-3, "color", "green", "marker", "s",
       phi4, [s4.My] * 1e-3, "color", "red", "marker", "s");
  legend("M_x", "M_y", "M_x", "M_y", "location", "northeast");
  grid;
  xlabel('\phi [°]'); ylabel('[Nm/mm]');

  print("M.svg", "-dsvg");

  fclose(fid);
