function [ point_coord, face_order, face_point ] = ...
  truncated_octahedron_shape_3d ( point_num, face_num, face_order_max )

%*****************************************************************************80
%
%% truncated_octahedron_shape_3d() describes a truncated octahedron in 3D.
%
%  Discussion:
%
%    The shape is a truncated octahedron.  There are 8 hexagons and 6
%    squares.
%
%    The truncated octahedron is an interesting shape because it
%    is "space filling".  In other words, all of 3D space can be
%    filled by a regular lattice of these shapes.
%
%    Call TRUNCATED_OCTAHEDRON_SIZE_3D to get the values of POINT_NUM,
%    FACE_NUM, and FACE_ORDER_MAX, so you can allocate space for the arrays.
%
%    For each face, the face list must be of length FACE_ORDER_MAX.
%    In cases where a face is of lower than maximum order (the
%    squares, in this case), the extra entries are listed as "-1".
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    15 August 2005
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer POINT_NUM, the number of points in the shape (24).
%
%    integer FACE_NUM, the number of faces in the shape (14).
%
%    integer FACE_ORDER_MAX, the maximum order of any face (6).
%
%  Output:
%
%    real POINT_COORD(3,POINT_NUM), the vertices.
%
%    integer FACE_ORDER(FACE_NUM), the number of vertices per face.
%
%    integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM); FACE_POINT(I,J)
%    contains the index of the I-th point in the J-th face.  The
%    points are listed in the counter-clockwise direction defined
%    by the outward normal at the face.
%
  dim_num = 3;
%
%  Set the point coordinates.
%
  point_coord(1:dim_num,1:point_num) = [  ...
    -1.5, -0.5,  0.0;        ...
    -1.5,  0.5,  0.0;        ...
    -1.0, -1.0, -0.70710677; ...
    -1.0, -1.0,  0.70710677; ...
    -1.0,  1.0, -0.70710677; ...
    -1.0,  1.0,  0.70710677; ...
    -0.5, -1.5,  0.0;        ...
    -0.5, -0.5, -1.4142135;  ...
    -0.5, -0.5,  1.4142135;  ...
    -0.5,  0.5, -1.4142135;  ...
    -0.5,  0.5,  1.4142135;  ...
    -0.5,  1.5,  0.0;        ...
     0.5, -1.5,  0.0;        ...
     0.5, -0.5, -1.4142135;  ...
     0.5, -0.5,  1.4142135;  ...
     0.5,  0.5, -1.4142135;  ...
     0.5,  0.5,  1.4142135;  ...
     0.5,  1.5,  0.0;        ...
     1.0, -1.0, -0.70710677; ...
     1.0, -1.0,  0.70710677; ...
     1.0,  1.0, -0.70710677; ...
     1.0,  1.0,  0.70710677; ... 
     1.5, -0.5,  0.0;        ...
     1.5,  0.5,  0.0 ]';
%
%  Set the face orders.
%
  face_order(1:face_num) = [ 4, 4, 4, 4, 4, 4, 6, 6, 6, 6, 6, 6, 6, 6 ];
%
%  Set faces.
%
  face_point(1:face_order_max,1:face_num) = [ ...
    17, 11,  9, 15, -1, -1; ...
    14,  8, 10, 16, -1, -1; ...
    22, 24, 21, 18, -1, -1; ...
    12,  5,  2,  6, -1, -1; ...
    13, 19, 23, 20, -1, -1; ...
     4,  1,  3,  7, -1, -1; ...
    19, 13,  7,  3,  8, 14; ...
    15,  9,  4,  7, 13, 20; ...
    16, 10,  5, 12, 18, 21; ...
    22, 18, 12,  6, 11, 17; ...
    20, 23, 24, 22, 17, 15; ...
    14, 16, 21, 24, 23, 19; ...
     9, 11,  6,  2,  1,  4; ...
     3,  1,  2,  5, 10,  8 ]';

  return
end
