Geometry/cube_shape_3d.f90 at master · lguez/Geometry · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
subroutine cube_shape_3d ( point_num, face_num, face_order_max, &
  point_coord, face_order, face_point )

!*****************************************************************************80
!
!! CUBE_SHAPE_3D describes a cube in 3D.
!
!  Discussion:
!
!    The vertices lie on the unit sphere.
!
!    The dual of the cube is the octahedron.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    05 October 2003
!
!  Author:
!
!    John Burkardt
!
!  Parameters:
!
!    Input, integer POINT_NUM, the number of points.
!
!    Input, integer FACE_NUM, the number of faces.
!
!    Input, integer FACE_ORDER_MAX, the maximum number of vertices
!    in a face.
!
!    Output, double precision POINT_COORD(3,POINT_NUM),
!    the vertices.
!
!    Output, integer FACE_ORDER(FACE_NUM), the number of vertices
!    per face.
!
!    Output, 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.
!
  implicit none

  integer face_num
  integer face_order_max
  integer, parameter :: dim_num = 3
  integer point_num

  double precision a
  integer face_order(face_num)
  integer face_point(face_order_max,face_num)
  double precision point_coord(dim_num,point_num)
!
!  Set point coordinates.
!
  a = sqrt ( 1.0D+00 / 3.0D+00 )

  point_coord(1:dim_num,1:point_num) = reshape ( (/ &
     -a, -a, -a, &
      a, -a, -a, &
      a,  a, -a, &
     -a,  a, -a, &
     -a, -a,  a, &
      a, -a,  a, &
      a,  a,  a, &
     -a,  a,  a /), (/ dim_num, point_num /) )
!
!  Set the face orders.
!
  face_order(1:face_num) = (/ &
    4, 4, 4, 4, 4, 4 /)
!
!  Set the faces.
!
  face_point(1:face_order_max,1:face_num) = reshape ( (/ &
     1, 4, 3, 2, &
     1, 2, 6, 5, &
     2, 3, 7, 6, &
     3, 4, 8, 7, &
     1, 5, 8, 4, &
     5, 6, 7, 8 /), (/ face_order_max, face_num /) )

  return
end