# Example: 2-dimensional truss
#
# -------------------------------------------------------------------

# Define output device depending on operating system

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

# Data

  F     =   1000;  % Force in N
  alpha =     30;  % Direction of force in degrees
  E     = 210000;  % Young's modulus in MPa
  A     =    400;  % Cross section area in mm^2
  a     =   1000;  % Length in mm

# Output file

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

# Model definition
# ----------------

# Model type and subtype

  model = struct("type", "solid", "subtype", "2d");

# Nodal points

  nodes( 1).id =  1; nodes( 1).coor = [ 0, 0];
  nodes( 2).id =  2; nodes( 2).coor = [ a, 0];
  nodes( 3).id =  3; nodes( 3).coor = [ 0, a];
  nodes( 4).id =  4; nodes( 4).coor = [ a, a];
  nodes( 5).id =  5; nodes( 5).coor = [ 0, 2 * a];
  nodes( 6).id =  6; nodes( 6).coor = [ a, 2 * a];
  nodes( 7).id =  7; nodes( 7).coor = [ 2 * a, 2 * a];
  nodes( 8).id =  8; nodes( 8).coor = [ 3 * a, 2 * a];
  nodes( 9).id =  9; nodes( 9).coor = [ 0.5 * a, 2.5 * a];
  nodes(10).id = 10; nodes(10).coor = [ 1.5 * a, 2.5 * a];
  nodes(11).id = 11; nodes(11).coor = [ 2.5 * a, 2.5 * a];

  model.nodes = nodes;

# Geometry and material

  geom = struct("A", A);
  mat  = struct("type", "iso", "E", E);

# Elements

  rod( 1).id = 1; rod( 1).nodes = [  1,  3];
  rod( 2).id = 2; rod( 2).nodes = [  1,  4];
  rod( 3).id = 3; rod( 3).nodes = [  2,  4];
  rod( 4).id = 4; rod( 4).nodes = [  3,  4];
  rod( 5).id = 5; rod( 5).nodes = [  3,  5];
  rod( 6).id = 6; rod( 6).nodes = [  4,  5];
  rod( 7).id = 7; rod( 7).nodes = [  4,  6];
  rod( 8).id = 8; rod( 8).nodes = [  5,  6];
  rod( 9).id = 9; rod( 9).nodes = [  6,  7];
  rod(10).id =10; rod(10).nodes = [  7,  8];
  rod(11).id =11; rod(11).nodes = [  5,  9];
  rod(12).id =12; rod(12).nodes = [  9,  6];
  rod(13).id =13; rod(13).nodes = [  6, 10];
  rod(14).id =14; rod(14).nodes = [ 10,  7];
  rod(15).id =15; rod(15).nodes = [  7, 11];
  rod(16).id =16; rod(16).nodes = [ 11,  8];
  rod(17).id =17; rod(17).nodes = [  9, 10];
  rod(18).id =18; rod(18).nodes = [ 10, 11];

  for n = 1 : length(rod)
      rod(n).type = "r2";
      rod(n).geom = geom; rod(n).mat  = mat;
  endfor

  model.elements = rod;

# Constraints

  hinge(1) = struct("id", 1, "dofs", [1, 2]);
  hinge(2) = struct("id", 2, "dofs", [1, 2]);

  model.constraints.prescribed = hinge;

# Load

  force = struct("id", 8,
                 "data", F * [sind(alpha), -cosd(alpha)]);

  model.loads.point = force;

# Analysis
# --------

# Create and plot component

  truss = mfs_new(fid, model);
  mfs_plot(truss, "nodid", 1, "eltid", 1, "force", 2,
           "fontsize", 18, "paperposition", [0, 0, 14, 10]);
  print(["model", EXT] , FORMAT, "-F:10");
  mfs_print(fid, truss, "loads", "point");

# Compute displacements

  truss = mfs_stiff(truss);
  truss = mfs_statresp(truss);
  ES = mfs_getresp(truss, "statresp", "strnegy");
  WR = mfs_getresp(truss, "statresp", "workrsl");
  fprintf(fid, "Strain energy   : %10.3e\n", ES);
  fprintf(fid, "Work of residual: %10.3e\n", WR);

  mfs_print(fid, truss, "statresp", "disp", "reac");
  mfs_plot(truss, "deform", 1, "nodid", 1, "force", 2,
           "reac", 2, "fontsize", 18);
  print(["freebody", EXT], FORMAT, "-F:10");

# Compute element results

  truss = mfs_results(truss, "statresp", "element");
  mfs_print(fid, truss, "statresp",
            "resultant", "stress", "strain");

# Export results to Gmsh

  mfs_export("truss_a.msh", "msh", truss, "mesh");
  mfs_export("truss_a.pos", "msh", truss,
             "statresp", "resultant", "stress");

  fclose(fid);
