Neko  0.9.99
A portable framework for high-order spectral element flow simulations
curve.f90
Go to the documentation of this file.
1 
2 module curve
3  use num_types, only : dp
4  use structs, only :struct_curve_t
5  use stack, only : stack_curve_t
6  use utils
7 
8  implicit none
9  private
10 
11  ! Maybe should be moved somewhere else
12  type, public :: curve_t
13  type(struct_curve_t), allocatable :: curve_el(:)
14  type(stack_curve_t), private :: scratch
15  integer :: size = 0
16  logical, private :: finalized = .false.
17  contains
18  procedure, pass(z) :: init => curve_element_init
19  procedure, pass(z) :: free => curve_element_free
20  procedure, pass(z) :: finalize => curve_element_finalize
21  procedure, pass(z) :: add_element => curve_element_add
22 ! procedure, pass(z) :: apply_xyz => curve_apply_xyz
23  end type curve_t
24 
25 contains
26 
28  subroutine curve_element_init(z, size)
29  class(curve_t), intent(inout) :: z
30  integer, optional :: size
31 
32  call curve_element_free(z)
33 
34  if (present(size)) then
35  call z%scratch%init(size)
36  else
37  call z%scratch%init()
38  end if
39 
40  end subroutine curve_element_init
41 
43  subroutine curve_element_free(z)
44  class(curve_t), intent(inout) :: z
45  if (allocated(z%curve_el)) then
46  deallocate(z%curve_el)
47  end if
48 
49  z%finalized = .false.
50  z%size = 0
51 
52  call z%scratch%free()
53 
54  end subroutine curve_element_free
55 
58  subroutine curve_element_finalize(z)
59  class(curve_t), target, intent(inout) :: z
60  integer :: i
61 
62  if (.not. z%finalized) then
63 
64  allocate(z%curve_el(z%scratch%size()))
65  select type (tp=>z%scratch%data)
66  type is (struct_curve_t)
67  do i = 1, z%scratch%size()
68  z%curve_el(i) = tp(i)
69  end do
70  end select
71 
72  z%size = z%scratch%size()
73 
74  call z%scratch%clear()
75 
76  z%finalized = .true.
77 
78  end if
79 
80  end subroutine curve_element_finalize
81 
83  subroutine curve_element_add(z, el_idx, curve_data, curve_type )
84  class(curve_t), intent(inout) :: z
85  real(kind=dp), dimension(5,12), intent(in) :: curve_data
86  integer, dimension(12), intent(in) :: curve_type
87  integer, intent(in) :: el_idx
88  type(struct_curve_t) :: c_el
89 
90  if (z%finalized) then
91  call neko_error('Zone already finalized')
92  end if
93  c_el%curve_data = curve_data
94  c_el%curve_type = curve_type
95  c_el%el_idx = el_idx
96  call z%scratch%push(c_el)
97  end subroutine curve_element_add
98 
99 ! !> Add a (facet, el) tuple to an unfinalized domain
100 ! subroutine curve_apply_xyz(this, x, y, z, Xh, elements, nel)
101 ! class(curve_t), intent(inout) :: this
102 ! integer, intent(in) :: lx, nel
103 ! class(element_t)
104 ! integer :: el_idx
105 ! real(kind=dp), dimension(Xh%lxyz, nel) :: x, y, z
106 ! do i = 1, this%size
107 ! el_idx = this%curve_el(i)%el_idx
108 ! do j = 1, 12
109 ! if (this%curve_el(i)%curve_type(j) .eq. 1) then
110 ! !call sphere_surface(j, this%curve_el(i)%curve_data(:,j),x(1,el_idx), y(1,el_idx), z(1, el_idx))
111 ! else if (this%curve_el(i)%curve_type(j) .eq. 2) then
112 ! !call generate_surface(j, this%curve_el(i)%curve_data(:,j), x(1,1,1,el_idx), y(1,1,1,el_idx), z(1,1,1, el_idx))
113 ! end if
114 ! end do
115 ! do j = 1, 12
116 ! if (this%curve_el(i)%curve_type(j) .eq. 3) then
117 ! call arc_surface(j, this%curve_el(i)%curve_data(:,j),x(1,1,1,el_idx), y(1,1,1,el_idx), z(1,1,1, el_idx))
118 ! end if
119 ! end do
120 !
121 ! end do
122 ! end subroutine curve_apply_xyz
123 
124 
125 ! subroutine sphere_surface(face, curve_data, x, y, z, Xh)
126 !! Major overhaul.
127 !! Martin Karp 1 March 2021 15:30:15
128 !!
129 !! 5 Aug 1988 19:29:52
130 !!
131 !! Program to generate spherical shell elements for NEKTON
132 !! input. Paul F. Fischer
133 !
134 ! NXM1 = NX-1
135 ! NYM1 = NY-1
136 ! NXY = NX*NZ
137 ! NXY3 = 3*NX*NZ
138 ! XCTR = curve_data(1,face)
139 ! YCTR = curve_data(2,face)
140 ! ZCTR = curve_data(3,face)
141 ! RADIUS = curve_data(4,face)
142 ! IFACE = EFACE1(IFCE)
143 !
144 ! Generate (normalized) corner vectors XCV(1,i,j):
145 !
146 ! CALL CRN3D(XCV,XC(1,IE),YC(1,IE),ZC(1,IE),CURVE(1,IFCE,IE),IFACE)
147 !
148 ! Generate edge vectors on the sphere RR=1.0,
149 ! for (r,s) = (-1,*),(1,*),(*,-1),(*,1)
150 !
151 ! CALL EDG3D(XYSRF,XCV(1,1,1),XCV(1,1,2), 1, 1, 1,NY,NX,NY)
152 ! CALL EDG3D(XYSRF,XCV(1,2,1),XCV(1,2,2),NX,NX, 1,NY,NX,NY)
153 ! CALL EDG3D(XYSRF,XCV(1,1,1),XCV(1,2,1), 1,NX, 1, 1,NX,NY)
154 ! CALL EDG3D(XYSRF,XCV(1,1,2),XCV(1,2,2), 1,NX,NY,NY,NX,NY)
155 !
156 ! ivtx = cface(1,ifce)
157 ! ivto = cface(2,ifce)
158 ! vout(1) = xc(ivtx,ie)-xc(ivto,ie)
159 ! vout(2) = yc(ivtx,ie)-yc(ivto,ie)
160 ! vout(3) = zc(ivtx,ie)-zc(ivto,ie)
161 !
162 ! vsph(1) = xc(ivtx,ie)-xctr
163 ! vsph(2) = yc(ivtx,ie)-yctr
164 ! vsph(3) = zc(ivtx,ie)-zctr
165 ! ifconcv = .true.
166 ! sign = DOT(vsph,vout,3)
167 ! if (sign.gt.0) ifconcv = .false.
168 ! write(6,*) 'THIS IS SIGN:',sign
169 !
170 ! DO 200 J=2,NYM1
171 ! CALL CROSS(VN1,XYSRF(1,1,J),XYSRF(1,NX,J))
172 ! DO 200 I=2,NXM1
173 ! CALL CROSS(VN2,XYSRF(1,I,1),XYSRF(1,I,NY))
174 ! if (ifconcv) then
175 ! IF (IFACE.EQ.1.OR.IFACE.EQ.4.OR.IFACE.EQ.5) THEN
176 ! CALL CROSS(XYSRF(1,I,J),VN2,VN1)
177 ! ELSE
178 ! CALL CROSS(XYSRF(1,I,J),VN1,VN2)
179 ! ENDIF
180 !20 CONTINUE
181 !
182 ! Normalize all vectors to the unit sphere.
183 !
184 ! DO 300 I=1,NXY
185 ! CALL NORM3D(XYSRF(1,I,1))
186 !30 CONTINUE
187 !
188 ! Scale by actual radius
189 !
190 ! CALL CMULT(XYSRF,RADIUS,NXY3)
191 !
192 ! Add back the sphere center offset
193 !
194 ! DO 400 I=1,NXY
195 ! XYSRF(1,I,1)=XYSRF(1,I,1)+XCTR
196 ! XYSRF(2,I,1)=XYSRF(2,I,1)+YCTR
197 ! XYSRF(3,I,1)=XYSRF(3,I,1)+ZCTR
198 !40 CONTINUE
199 !
200 !
201 ! Transpose data, if necessary
202 !
203 ! IF (IFACE.EQ.1.OR.IFACE.EQ.4.OR.IFACE.EQ.5) THEN
204 ! DO 500 J=1 ,NY
205 ! DO 500 I=J+1,NX
206 ! TMP=XYSRF(1,I,J)
207 ! XYSRF(1,I,J)=XYSRF(1,J,I)
208 ! XYSRF(1,J,I)=TMP
209 ! TMP=XYSRF(2,I,J)
210 ! XYSRF(2,I,J)=XYSRF(2,J,I)
211 ! XYSRF(2,J,I)=TMP
212 ! TMP=XYSRF(3,I,J)
213 ! XYSRF(3,I,J)=XYSRF(3,J,I)
214 ! XYSRF(3,J,I)=TMP
215 !50 CONTINUE
216 ! ENDIF
217 !
218 ! Compute surface deflection and perturbation due to face IFACE
219 !
220 ! CALL DSSET(NX,NY,NZ)
221 ! JS1 = SKPDAT(1,IFACE)
222 ! JF1 = SKPDAT(2,IFACE)
223 ! JSKIP1 = SKPDAT(3,IFACE)
224 ! JS2 = SKPDAT(4,IFACE)
225 ! JF2 = SKPDAT(5,IFACE)
226 ! JSKIP2 = SKPDAT(6,IFACE)
227 !
228 ! IOPP(1) = NX-1
229 ! IOPP(2) = NX*(NY-1)
230 ! IOPP(3) = NX*NY*(NZ-1)
231 ! NXX(1) = NX
232 ! NXX(2) = NY
233 ! NXX(3) = NZ
234 ! IDIR = 2*MOD(IFACE,2) - 1
235 ! IFC2 = (IFACE+1)/2
236 ! DELT = 0.0
237 ! I=0
238 ! DO 700 J2=JS2,JF2,JSKIP2
239 ! DO 700 J1=JS1,JF1,JSKIP1
240 ! I=I+1
241 ! JOPP = J1 + IOPP(IFC2)*IDIR
242 ! X2(1) = XML(J1,J2,1,IE)
243 ! X2(2) = YML(J1,J2,1,IE)
244 ! X2(3) = ZML(J1,J2,1,IE)
245 !
246 ! DX(1) = XYSRF(1,I,1)-X2(1)
247 ! DX(2) = XYSRF(2,I,1)-X2(2)
248 ! DX(3) = XYSRF(3,I,1)-X2(3)
249 !
250 ! NXS = NXX(IFC2)
251 ! JOFF = (J1-JOPP)/(NXS-1)
252 ! DO 600 IX = 2,NXS
253 ! J = JOPP + JOFF*(IX-1)
254 ! ZETA = 0.5*(ZGML(IX,IFC2) + 1.0)
255 ! XML(J,J2,1,IE) = XML(J,J2,1,IE)+DX(1)*ZETA
256 ! YML(J,J2,1,IE) = YML(J,J2,1,IE)+DX(2)*ZETA
257 ! ZML(J,J2,1,IE) = ZML(J,J2,1,IE)+DX(3)*ZETA
258 !60 CONTINUE
259 !70 CONTINUE
260 !
261 ! return
262 ! end
263 
264 
265 end module curve
Defines a domain as a subset of facets in a mesh.
Definition: curve.f90:2
subroutine curve_element_free(z)
Deallocate a domain.
Definition: curve.f90:44
subroutine curve_element_finalize(z)
Finalize a domain list.
Definition: curve.f90:59
subroutine curve_element_init(z, size)
Initialize a curved domain.
Definition: curve.f90:29
subroutine curve_element_add(z, el_idx, curve_data, curve_type)
Add a (facet, el) tuple to an unfinalized domain.
Definition: curve.f90:84
integer, parameter, public dp
Definition: num_types.f90:9
Implements a dynamic stack ADT.
Definition: stack.f90:35
Defines structs that are used... Dont know if we should keep it though.
Definition: structs.f90:2
Utilities.
Definition: utils.f90:35
Curved element stack.
Definition: stack.f90:112