# Example: Gull wing
#
# Solid model:
#  - Definition of the solid model
#  - Computation of first five normal modes to check the model
#  - Storage of solid component including node sets in binary file
#    solid.bin
#
# ---------------------------------------------------------------------

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

# Data (N, mm, s)
# ---------------

# Geometry

  x1  =    100;  % x-position of front spar
  x2  =    375;  % x-position of main spar
  x3  =   1200;  % x-position of rear spar at wing root
  x4  =    800;  % x-position of rear spar at wing tip

  ya  =   4000;  % y-position of kink
  yt  =   6000;  % y-position of wing tip
  zt  =   1000;  % z-position of wing tip

# Material

  E   =  70000;  % Young's modulus 
  ny  =   0.34;  % Poisson ratio
  rho = 2.7E-9;  % Mass density

# Cross sections

  frontspar = struct("type", "I",               % Front Spar
                     "b",     30, "h", 100,
                     "t",      2, "s",   2);

  mainspar = struct("type", "box",              % Main Spar
                    "styp", "thin", "b", 150,
                    "h",       180, "t",   2);

  rearspar = struct("type", "I",                % Rear Spar
                    "b",     30, "h", 100,
                    "t",      2, "s",   2);

  ribs     = struct("type", "I",                % Ribs
                    "b",     40, "h", 120,
                    "t",      1, "s",   1);
  
# Model Definition
# ----------------

# Model type

  solid = struct("type", "solid", "subtype", "3d");

# Material

  mat = struct("type", "iso", "E", E, "ny", ny, "rho", rho);

# Front spar

  geom   = mfs_beamsection(frontspar.type, frontspar.b, 
                           frontspar.h, frontspar.t, 
                           frontspar.s);
  geom.v = [0, 0, 1];

  nodes(1) = struct("id",  101, "coor", [ 100,    0,    0]);
  nodes(2) = struct("id",  901, "coor", [ 100, 4000,    0]);
  nodes(3) = struct("id", 1301, "coor", [ 100, 6000, 1000]);

  idnew = 201 : 100 : 801;
  idelt = 1 : 8;

  [nodes, elem1] = mfs_line(nodes, 101, 901, idnew, idelt, 
                            "b2", geom, mat);

  idnew = 1001 : 100 : 1201;
  idelt = 9 : 12;

  [nodes,elem2] = mfs_line(nodes, 901, 1301, idnew, idelt, 
                           "b2", geom, mat);

# Main spar

  geom   = mfs_beamsection(mainspar.type, mainspar.styp,
                           mainspar.b, mainspar.h,
                           mainspar.t);
  geom.v = [0, 0, 1];

  n = length(nodes);

  nodes(++n) = struct("id",  111, "coor", [ 375,    0,    0]);
  nodes(++n) = struct("id",  911, "coor", [ 375, 4000,    0]);
  nodes(++n) = struct("id", 1311, "coor", [ 375, 6000, 1000]);

  idnew = 211 : 100 : 811;
  idelt = 13 : 20;

  [nodes, elem3] = mfs_line(nodes, 111, 911, idnew, idelt, 
                            "b2", geom, mat);

  idnew = 1011 : 100 : 1211;
  idelt = 21 : 24;

  [nodes,elem4] = mfs_line(nodes, 911, 1311, idnew, idelt, 
                           "b2", geom, mat);

# Rear spar

  geom   = mfs_beamsection(rearspar.type, rearspar.b,
                           rearspar.h, rearspar.t,
                           rearspar.s);
  geom.v = [0, 0, 1];

  n = length(nodes);

  nodes(++n) = struct("id",  141, "coor", [ 1200,    0,    0]);
  nodes(++n) = struct("id",  941, "coor", [ 1200, 4000,    0]);
  nodes(++n) = struct("id", 1341, "coor", [  800, 6000, 1000]);

  idnew = 241 : 100 : 841;
  idelt = 25 : 32;

  [nodes, elem5] = mfs_line(nodes, 141, 941, idnew, idelt, 
                            "b2", geom, mat);

  idnew = 1041 : 100 : 1241;
  idelt = 33 : 36;

  [nodes,elem6] = mfs_line(nodes, 941, 1341, idnew, idelt, 
                           "b2", geom, mat);


  elemh = [elem1, elem2, elem3, elem4, elem5, elem6];

# Ribs in front of main spar

  geomi   = mfs_beamsection(ribs.type, ribs.b, ribs.h,
                            ribs.t, ribs.s);
  geomi.v = [0, 0, 1];

  geoma   = mfs_beamsection(ribs.type, ribs.b, ribs.h,
                            ribs.t, ribs.s);
  geoma.v = [0, -1 , 2];

  ide = length(elemh) + 1;
  id1 = 101; id2 = 111;

  for n = 1 : 9
      elemr(n) = struct("id",   ide++, "nodes", [id1, id2],
                        "type",  "b2", "geom",  geomi,
                        "mat",    mat);
      id1 += 100; id2 += 100;
  endfor

  for n = 10 : 13
      elemr(n) = struct("id",   ide++, "nodes", [id1, id2],
                        "type",  "b2", "geom",  geoma,
                        "mat",    mat);
      id1 += 100; id2 += 100;
  endfor

# Ribs behind main spar

  id1 = 111; id2 = 141; idnew = [121, 131]; 

  for n = 1 : 9
      idelt = ide : ide + 2;
      [nodes, elemx] = mfs_line(nodes, id1, id2, idnew, idelt,
                                 "b2", geomi, mat);
      id1 += 100; id2 += 100; idnew += 100;
      ide += 3;
      elemr = [elemr, elemx];
  endfor

  for n = 1 : 4
      idelt = ide : ide + 2;
      [nodes, elemx] = mfs_line(nodes, id1, id2, idnew, idelt,
                                 "b2", geoma, mat);
      id1 += 100; id2 += 100; idnew += 100;
      ide += 3;
      elemr = [elemr, elemx];
  endfor

  solid.nodes     = nodes;
  solid.elements  = [elemh, elemr];

# Node sets (needed to define the splines)

  nset.inner_wing = [ 101 : 100 : 901, 111 : 100 : 911, ...
                      121 : 100 : 911, 131 : 100 : 931, ...
                      141 : 100 : 941 ];
  nset.outer_wing = [ 901 : 100 : 1301, 911 : 100 : 1311, ...
                      921 : 100 : 1321, 931 : 100 : 1331, ...
                      941 : 100 : 1341 ];

  solid.nset = nset;

# Constraints

  prescribed = struct("id",   {  101,   111,   141},
                      "dofs", {1 : 3, 1 : 6, 1 : 3});
  solid.constraints.prescribed = prescribed;

# Analysis
# --------

# Create and export component

  wings = mfs_new(fid, solid);
  mfs_export("solid.msh", "msh", wings, "mesh", "mesh", "axes");

# Compute stiffness and mass matrices

  wings = mfs_stiff(wings);
  wings = mfs_mass(wings);
  mfs_massproperties(fid, wings);

# Save the solid component

  save solid.bin wings

# Compute and export the normal modes

  wings = mfs_freevib(wings, 5);
  mfs_print(fid, wings, "modes", "freq");
  mfs_export("modes.dsp", "msh", wings, "modes", "disp");

  fclose(fid);
