function mdl = model(ne, etype)

# Input  ne        Number of elements along one direction
#        etype     Element name
# Output mdl       Structure with model definition
#
# The function generates the model of the Scordelis-Lo roof.
#
# Node and element numbers increase by one in circumferential
# direction and by row increments in longitudinal direction.
#
# ---------------------------------------------------------------

# Data (N, mm)

  R    =     7620;  % Radius
  L    =     7620;  % Length of model
  t    =     76.2;  % Thickness
  phiE =       40;  % Angle in degress

  E    = 2.0684e4;  % Young's modulus
  ny   =        0;  % Poisson's ratio

  pz   = 4.309e-3;  % Weight per area
  g    =     9810;  % Gravity acceleration

# Nodal point data

  np   = ne + 1;    % Nodes in one row
  incp = np;        % Increment of nodal point identifiers

# Model type and subtype

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

# Geometry and material

  geom = struct("t", t);
  mat  = struct("type", "iso", "E", E, "ny", ny,
                "rho", pz / (t * g));

# Nodes

  phi = linspace(0, phiE * pi / 180, np);
  x   = R * sin(phi);
  y   = linspace(0, L, np);
  z   = R * cos(phi);

  k = 1;

  for n = 1 : np
      for m = 1 : np
          nodes(k) = struct("id", k++,
                            "coor", [x(m), y(n), z(m)]);
      end
  end

  mdl.nodes = nodes;

  nofnod = np * np;
  mdl.nset.FG = (nofnod - np + 1) : nofnod; % Nodes along FG
  mdl.nset.G = mdl.nset.FG(end);            % Point G

# Elements

  k = 1;

  switch etype

   case "s4"

     enode1 = [1, 2, 2 + incp, 1 + incp];

     for n = 1 : ne
         enodes = enode1;
         for m = 1 : ne
             elem(k) = struct("id", k++, "type", etype,
                              "nodes", enodes++,
                              "geom", geom, "mat", mat);
         endfor
         enode1 += incp;
     endfor

     mdl.eset.FG = (elem(end).id - ne + 1) : elem(end).id;

   case "s3"

     enode1 = [1, 2, 1 + incp]; enode2 = [2 + incp, 1 + incp, 2];

     for n = 1 : ne
         enodea = enode1; enodeb = enode2;
         for m = 1 : ne
             elem(k) = struct("id", k++, "type", etype,
                              "nodes", enodea++,
                              "geom", geom, "mat", mat);
             elem(k) = struct("id", k++, "type", etype,
                              "nodes", enodeb++,
                              "geom", geom, "mat", mat);
         endfor
         enode1 += incp; enode2 += incp;
     endfor

     FG = (elem(end).id - 2 * ne + 1) : 2 : elem(end).id;
     mdl.eset.FG = FG;

  endswitch

  mdl.elements = elem;

# Constraints

  k = 1;

  % at y = 0: u = w = 0

  for m = 1 : np
      fix(k++) = struct("id", m, "dofs", [1, 3]);
  endfor

  % at phi = 0: symmetry

  n2 = 1 + ne * incp;
  for n = 1 : incp : n2
      fix(k++) = struct("id", n, "dofs", [1, 5, 6]);
  endfor

  % at y = L: symmetry

  m1 = n2; m2 = m1 + ne;
  for m = m1 : m2
      fix(k++) = struct("id", m, "dofs", [2, 4, 6]);
  endfor

  mdl.constraints.prescribed = fix;

# Load

  mdl.loads.inertia = struct("data", [0, 0, g]);

endfunction
