Skip to content

Commit 0cabfff

Browse files
danieljvickersdanieljvickers
committed
Finished adding support for existing IB geometries. They will all now shift and rotate the grid cell into their local frame and then compute validity from the local frame. This is true for all IB patches except for circles and spheres, as rotational symetry solves the problem for me in those cases.
1 parent ddfcd0f commit 0cabfff

File tree

1 file changed

+78
-109
lines changed

1 file changed

+78
-109
lines changed

‎src/common/m_ib_patches.fpp‎

Lines changed: 78 additions & 109 deletions
Original file line numberDiff line numberDiff line change
@@ -180,7 +180,7 @@ contains
180180
integer, intent(in) :: patch_id
181181
integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
182182

183-
real(wp) :: x0, y0, f, x_act, y_act, ca_in, pa, ma, ta, theta
183+
real(wp) :: f, ca_in, pa, ma, ta, theta
184184
real(wp) :: xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
185185
integer :: i, j, k
186186
integer :: Np1, Np2
@@ -207,26 +207,26 @@ contains
207207
#endif
208208
Np = Np1 + Np2 + 1
209209

210-
allocate (airfoil_grid_u(1:Np))
210+
allocate (airfoil_grid_u(1:Np))
211211
allocate (airfoil_grid_l(1:Np))
212212

213-
airfoil_grid_u(1)%x = x0
214-
airfoil_grid_u(1)%y = y0
213+
airfoil_grid_u(1)%x = 0._wp
214+
airfoil_grid_u(1)%y = 0._wp
215215

216-
airfoil_grid_l(1)%x = x0
217-
airfoil_grid_l(1)%y = y0
216+
airfoil_grid_l(1)%x = 0._wp
217+
airfoil_grid_l(1)%y = 0._wp
218218

219219
eta = 1._wp
220220

221221
do i = 1, Np1 + Np2 - 1
222222
if (i <= Np1) then
223-
xc = x0 + i*(pa*ca_in/Np1)
224-
xa = (xc - x0)/ca_in
223+
xc = i*(pa*ca_in/Np1)
224+
xa = xc/ca_in
225225
yc = (ma/pa**2)*(2*pa*xa - xa**2)
226226
dycdxc = (2*ma/pa**2)*(pa - xa)
227227
else
228-
xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2)
229-
xa = (xc - x0)/ca_in
228+
xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2)
229+
xa = xc/ca_in
230230
yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
231231
dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
232232
end if
@@ -241,11 +241,11 @@ contains
241241
xl = xa + yt*sin_c
242242
yl = yc - yt*cos_c
243243

244-
xu = xu*ca_in + x0
245-
yu = yu*ca_in + y0
244+
xu = xu*ca_in
245+
yu = yu*ca_in
246246

247-
xl = xl*ca_in + x0
248-
yl = yl*ca_in + y0
247+
xl = xl*ca_in
248+
yl = yl*ca_in
249249

250250
airfoil_grid_u(i + 1)%x = xu
251251
airfoil_grid_u(i + 1)%y = yu
@@ -255,65 +255,57 @@ contains
255255

256256
end do
257257

258-
airfoil_grid_u(Np)%x = x0 + ca_in
259-
airfoil_grid_u(Np)%y = y0
258+
airfoil_grid_u(Np)%x = ca_in
259+
airfoil_grid_u(Np)%y = 0._wp
260260

261-
airfoil_grid_l(Np)%x = x0 + ca_in
262-
airfoil_grid_l(Np)%y = y0
261+
airfoil_grid_l(Np)%x = ca_in
262+
airfoil_grid_l(Np)%y = 0._wp
263263

264264
do j = 0, n
265265
do i = 0, m
266-
if (.not. f_is_default(patch_ib(patch_id)%theta)) then
267-
x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta)
268-
y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta)
269-
else
270-
x_act = x_cc(i) -
271-
y_act = y_cc(j)
272-
end if
273-
274-
xyz_local = [x_act, y_act, 0._wp] ! get coordinate frame centered on IB
275-
xyz_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordiantes
266+
xy_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, 0._wp] ! get coordinate frame centered on IB
267+
xy_local = matmul(inverse_rotation, xy_local) ! rotate the frame into the IB's coordiantes
276268
277-
if (x_act >= x0 .and. x_act <= x0 + ca_in) then
278-
xa = (x_act - x0)/ca_in
269+
if (xy_local(1) >= 0._wp .and. xy_local(1) <= ca_in) then
270+
xa = xy_local(1)/ca_in
279271
if (xa <= pa) then
280272
yc = (ma/pa**2)*(2*pa*xa - xa**2)
281273
dycdxc = (2*ma/pa**2)*(pa - xa)
282274
else
283275
yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
284276
dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
285277
end if
286-
if (y_act >= y0) then
278+
if (xy_local(1) >= 0._wp) then
287279
k = 1
288-
do while (airfoil_grid_u(k)%x < x_act .and. k <= Np)
280+
do while (airfoil_grid_u(k)%x < xy_local(1) .and. k <= Np)
289281
k = k + 1
290282
end do
291-
if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then
292-
if (y_act <= airfoil_grid_u(k)%y) then
283+
if (f_approx_equal(airfoil_grid_u(k)%x, xy_local(1))) then
284+
if (xy_local(1) <= airfoil_grid_u(k)%y) then
293285
!!IB
294286
ib_markers_sf(i, j, 0) = patch_id
295287
end if
296288
else
297-
f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
298-
if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
289+
f = (airfoil_grid_u(k)%x - xy_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
290+
if (xy_local(1) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
299291
!!IB
300292
ib_markers_sf(i, j, 0) = patch_id
301293
end if
302294
end if
303295
else
304296
k = 1
305-
do while (airfoil_grid_l(k)%x < x_act)
297+
do while (airfoil_grid_l(k)%x < xy_local(1))
306298
k = k + 1
307299
end do
308-
if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then
309-
if (y_act >= airfoil_grid_l(k)%y) then
300+
if (f_approx_equal(airfoil_grid_l(k)%x, xy_local(1))) then
301+
if (xy_local(1) >= airfoil_grid_l(k)%y) then
310302
!!IB
311303
ib_markers_sf(i, j, 0) = patch_id
312304
end if
313305
else
314-
f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
306+
f = (airfoil_grid_l(k)%x - xy_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
315307
316-
if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
308+
if (xy_local(1) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
317309
!!IB
318310
ib_markers_sf(i, j, 0) = patch_id
319311
end if
@@ -323,16 +315,6 @@ contains
323315
end do
324316
end do
325317
326-
if (.not. f_is_default(patch_ib(patch_id)%theta)) then
327-
do i = 1, Np
328-
airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0
329-
airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0
330-
331-
airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0
332-
airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0
333-
end do
334-
end if
335-
336318
end subroutine s_ib_airfoil
337319
338320
!! @param patch_id is the patch identifier
@@ -343,19 +325,21 @@ contains
343325
integer, intent(in) :: patch_id
344326
integer, dimension(0:m, 0:n, 0:p), intent(inout) :: ib_markers_sf
345327
346-
real(wp) :: x0, y0, z0, lz, z_max, z_min, f, x_act, y_act, ca_in, pa, ma, ta, theta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
328+
real(wp) :: lz, z_max, z_min, f, ca_in, pa, ma, ta, xa, yt, xu, yu, xl, yl, xc, yc, dycdxc, sin_c, cos_c
347329
integer :: i, j, k, l
348330
integer :: Np1, Np2
349331
350-
x0 = patch_ib(patch_id)%x_centroid
351-
y0 = patch_ib(patch_id)%y_centroid
352-
z0 = patch_ib(patch_id)%z_centroid
332+
real(wp), dimension(1:3) :: xyz_local !< x, y, z coordinates in local IB frame
333+
real(wp), dimension(1:3, 1:3) :: inverse_rotation
334+
335+
x_centroid = patch_ib(patch_id)%x_centroid
336+
y_centroid = patch_ib(patch_id)%y_centroid
337+
z_centroid = patch_ib(patch_id)%z_centroid
353338
lz = patch_ib(patch_id)%length_z
354339
ca_in = patch_ib(patch_id)%c
355340
pa = patch_ib(patch_id)%p
356341
ma = patch_ib(patch_id)%m
357342
ta = patch_ib(patch_id)%t
358-
theta = pi*patch_ib(patch_id)%theta/180._wp
359343
360344
! rank(dx) is not consistent between pre_process and simulation. This IFDEF prevents compilation errors
361345
#ifdef MFC_PRE_PROCESS
@@ -370,26 +354,26 @@ contains
370354
allocate (airfoil_grid_u(1:Np))
371355
allocate (airfoil_grid_l(1:Np))
372356
373-
airfoil_grid_u(1)%x = x0
374-
airfoil_grid_u(1)%y = y0
357+
airfoil_grid_u(1)%x = 0._wp
358+
airfoil_grid_u(1)%y = 0._wp
375359
376-
airfoil_grid_l(1)%x = x0
377-
airfoil_grid_l(1)%y = y0
360+
airfoil_grid_l(1)%x = 0._wp
361+
airfoil_grid_l(1)%y = 0._wp
378362
379-
z_max = z0 + lz/2
380-
z_min = z0 - lz/2
363+
z_max = lz/2
364+
z_min = -lz/2
381365
382366
eta = 1._wp
383367
384368
do i = 1, Np1 + Np2 - 1
385369
if (i <= Np1) then
386-
xc = x0 + i*(pa*ca_in/Np1)
387-
xa = (xc - x0)/ca_in
370+
xc = i*(pa*ca_in/Np1)
371+
xa = xc/ca_in
388372
yc = (ma/pa**2)*(2*pa*xa - xa**2)
389373
dycdxc = (2*ma/pa**2)*(pa - xa)
390374
else
391-
xc = x0 + pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2)
392-
xa = (xc - x0)/ca_in
375+
xc = pa*ca_in + (i - Np1)*((ca_in - pa*ca_in)/Np2)
376+
xa = xc/ca_in
393377
yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
394378
dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
395379
end if
@@ -404,11 +388,11 @@ contains
404388
xl = xa + yt*sin_c
405389
yl = yc - yt*cos_c
406390
407-
xu = xu*ca_in + x0
408-
yu = yu*ca_in + y0
391+
xu = xu*ca_in
392+
yu = yu*ca_in
409393
410-
xl = xl*ca_in + x0
411-
yl = yl*ca_in + y0
394+
xl = xl*ca_in
395+
yl = yl*ca_in
412396
413397
airfoil_grid_u(i + 1)%x = xu
414398
airfoil_grid_u(i + 1)%y = yu
@@ -418,85 +402,70 @@ contains
418402
419403
end do
420404
421-
airfoil_grid_u(Np)%x = x0 + ca_in
422-
airfoil_grid_u(Np)%y = y0
405+
airfoil_grid_u(Np)%x = ca_in
406+
airfoil_grid_u(Np)%y = 0._wp
423407
424-
airfoil_grid_l(Np)%x = x0 + ca_in
425-
airfoil_grid_l(Np)%y = y0
408+
airfoil_grid_l(Np)%x = ca_in
409+
airfoil_grid_l(Np)%y = 0._wp
426410
427411
do l = 0, p
428-
if (z_cc(l) >= z_min .and. z_cc(l) <= z_max) then
429-
do j = 0, n
430-
do i = 0, m
412+
do j = 0, n
413+
do i = 0, m
414+
xyz_local = [x_cc(i) - x_centroid, y_cc(j) - y_centroid, z_cc(l) - z_centroid] ! get coordinate frame centered on IB
415+
xyz_local = matmul(inverse_rotation, xyz_local) ! rotate the frame into the IB's coordiantes
431416

432-
if (.not. f_is_default(patch_ib(patch_id)%theta)) then
433-
x_act = (x_cc(i) - x0)*cos(theta) - (y_cc(j) - y0)*sin(theta) + x0
434-
y_act = (x_cc(i) - x0)*sin(theta) + (y_cc(j) - y0)*cos(theta) + y0
435-
else
436-
x_act = x_cc(i)
437-
y_act = y_cc(j)
438-
end if
417+
if (xyz_local(3) >= z_min .and. xyz_local(3) <= z_max) then
439418

440-
if (x_act >= x0 .and. x_act <= x0 + ca_in) then
441-
xa = (x_act - x0)/ca_in
419+
if (xyz_local(1) >= 0._wp .and. xyz_local(1) <= ca_in) then
420+
xa = xyz_local(1)/ca_in
442421
if (xa <= pa) then
443422
yc = (ma/pa**2)*(2*pa*xa - xa**2)
444423
dycdxc = (2*ma/pa**2)*(pa - xa)
445424
else
446425
yc = (ma/(1 - pa)**2)*(1 - 2*pa + 2*pa*xa - xa**2)
447426
dycdxc = (2*ma/(1 - pa)**2)*(pa - xa)
448427
end if
449-
if (y_act >= y0) then
428+
if (xyz_local(2) >= 0._wp) then
450429
k = 1
451-
do while (airfoil_grid_u(k)%x < x_act)
430+
do while (airfoil_grid_u(k)%x < xyz_local(1))
452431
k = k + 1
453432
end do
454-
if (f_approx_equal(airfoil_grid_u(k)%x, x_act)) then
455-
if (y_act <= airfoil_grid_u(k)%y) then
433+
if (f_approx_equal(airfoil_grid_u(k)%x, xyz_local(1))) then
434+
if (xyz_local(2) <= airfoil_grid_u(k)%y) then
456435
!!IB
457436
ib_markers_sf(i, j, l) = patch_id
458437
end if
459438
else
460-
f = (airfoil_grid_u(k)%x - x_act)/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
461-
if (y_act <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
439+
f = (airfoil_grid_u(k)%x - xyz_local(1))/(airfoil_grid_u(k)%x - airfoil_grid_u(k - 1)%x)
440+
if (xyz_local(2) <= ((1._wp - f)*airfoil_grid_u(k)%y + f*airfoil_grid_u(k - 1)%y)) then
462441
!!IB
463442
ib_markers_sf(i, j, l) = patch_id
464443
end if
465444
end if
466445
else
467446
k = 1
468-
do while (airfoil_grid_l(k)%x < x_act)
447+
do while (airfoil_grid_l(k)%x < xyz_local(1))
469448
k = k + 1
470449
end do
471-
if (f_approx_equal(airfoil_grid_l(k)%x, x_act)) then
472-
if (y_act >= airfoil_grid_l(k)%y) then
450+
if (f_approx_equal(airfoil_grid_l(k)%x, xyz_local(1))) then
451+
if (xyz_local(2) >= airfoil_grid_l(k)%y) then
473452
!!IB
474453
ib_markers_sf(i, j, l) = patch_id
475454
end if
476455
else
477-
f = (airfoil_grid_l(k)%x - x_act)/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
456+
f = (airfoil_grid_l(k)%x - xyz_local(1))/(airfoil_grid_l(k)%x - airfoil_grid_l(k - 1)%x)
478457

479-
if (y_act >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
458+
if (xyz_local(2) >= ((1._wp - f)*airfoil_grid_l(k)%y + f*airfoil_grid_l(k - 1)%y)) then
480459
!!IB
481460
ib_markers_sf(i, j, l) = patch_id
482461
end if
483462
end if
484463
end if
485464
end if
486-
end do
465+
end if
487466
end do
488-
end if
489-
end do
490-
491-
if (.not. f_is_default(patch_ib(patch_id)%theta)) then
492-
do i = 1, Np
493-
airfoil_grid_l(i)%x = (airfoil_grid_l(i)%x - x0)*cos(theta) + (airfoil_grid_l(i)%y - y0)*sin(theta) + x0
494-
airfoil_grid_l(i)%y = -1._wp*(airfoil_grid_l(i)%x - x0)*sin(theta) + (airfoil_grid_l(i)%y - y0)*cos(theta) + y0
495-
496-
airfoil_grid_u(i)%x = (airfoil_grid_u(i)%x - x0)*cos(theta) + (airfoil_grid_u(i)%y - y0)*sin(theta) + x0
497-
airfoil_grid_u(i)%y = -1._wp*(airfoil_grid_u(i)%x - x0)*sin(theta) + (airfoil_grid_u(i)%y - y0)*cos(theta) + y0
498467
end do
499-
end if
468+
end do
500469

501470
end subroutine s_ib_3D_airfoil
502471

0 commit comments

Comments
 (0)