Geometry/dual_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
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
subroutine dual_shape_3d ( point_num, face_num, face_order_max, &
  point_coord, face_order, face_point, point_num2, face_num2, &
  face_order_max2, point_coord2, face_order2, face_point2 )

!*****************************************************************************80
!
!! DUAL_SHAPE_3D constructs the dual of a shape in 3D.
!
!  Licensing:
!
!    This code is distributed under the GNU LGPL license.
!
!  Modified:
!
!    23 August 2005
!
!  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
!    per face.
!
!    Input, double precision POINT_COORD(3,POINT_NUM), the points.
!
!    Input, integer FACE_ORDER(FACE_NUM), the number of vertices
!    per face.
!
!    Input, integer FACE_POINT(FACE_ORDER_MAX,FACE_NUM);
!    FACE_POINT(I,J) is 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.
!
!    Input, integer POINT_NUM2, the number of points in the dual.
!
!    Input, integer FACE_NUM2, the number of faces in the dual.
!
!    Input, integer FACE_ORDER_MAX2, the maximum number of
!    vertices per face in the dual.
!
!    Output, double precision POINT_COORD2(3,POINT_NUM2), the point
!    coordinates of the dual.
!
!    Output, integer FACE_ORDER2(FACE_NUM2), the number of
!    vertices per face.
!
!    Output, integer FACE_POINT2(FACE_ORDER_MAX2,FACE_NUM2),
!    the vertices of each face in the dual.
!
  implicit none

  integer face_num
  integer face_num2
  integer face_order_max
  integer face_order_max2
  integer, parameter :: dim_num = 3
  integer point_num
  integer point_num2

  integer col
  integer face
  integer face_order(face_num)
  integer face_order2(face_num2)
  integer face_point(face_order_max,face_num)
  integer face_point2(face_order_max2,face_num2)
  integer i
  integer inext
  integer iprev
  integer istop
  integer j
  integer k
  double precision norm
  double precision p(dim_num)
  double precision point_coord(dim_num,point_num)
  double precision point_coord2(dim_num,point_num2)
  integer row
!
!  This computation should really compute the center of gravity
!  of the face, in the general case.
!
!  We'll also assume the vertices of the original and the dual
!  are to lie on the unit sphere, so we can normalize the
!  position vector of the vertex.
!
  do face = 1, face_num

    p(1:dim_num) = 0.0D+00

    do j = 1, face_order(face)
      k = face_point(j,face)
      p(1:dim_num) = p(1:dim_num) + point_coord(1:dim_num,k)
    end do

    norm = sqrt ( sum ( p(1:dim_num)**2 ) )

    point_coord2(1:dim_num,face) = p(1:dim_num) / norm

  end do
!
!  Now build the face in the dual associated with each node FACE.
!
  do face = 1, face_num2
!
!  Initialize the order.
!
    face_order2(face) = 0
!
!  Find the first occurrence of FACE in an edge of polyhedron.
!
    call i4col_find_item ( face_order_max, face_num, face_point, &
      face, row, col )

    if ( row <= 0 ) then
      write ( *, '(a)' ) ' '
      write ( *, '(a)' ) 'DUAL_SHAPE_3D - Fatal error!'
      write ( *, '(a,i8)' ) '  Could not find an edge using node ', face
      stop 1
    end if
!
!  Save the following node as ISTOP.
!  When we encounter ISTOP again, this will mark the end of our search.
!
    i = row + 1
    if ( face_order(col) < i ) then
      i = 1
    end if

    istop = face_point(i,col)
!
!  Save the previous node as INEXT.
!
    do

      i = row - 1
      if ( i < 1 ) then
        i = i + face_order(col)
      end if

      inext = face_point(i,col)

      face_order2(face) = face_order2(face) + 1

      face_point2(face_order2(face),face) = col
!
!  If INEXT =/= ISTOP, continue.
!
      if ( inext == istop ) then
        exit
      end if
!
!  Set IPREV:= INEXT.
!
      iprev = inext
!
!  Search for the occurrence of the edge FACE-IPREV.
!
      call i4col_find_pair_wrap ( face_order_max, face_num, face_point, &
        face, iprev, row, col )

      if ( row <= 0 ) then
        write ( *, '(a)' ) ' '
        write ( *, '(a)' ) 'DUAL_SHAPE_3D - Fatal error!'
        write ( *, '(a,i8)' ) '  No edge from node ', iprev
        write ( *, '(a,i8)' ) '  to node ', face
        stop 1
      end if

    end do

  end do

  return
end