13. Structured blocks and overset grids

The chapters so far put a body into a Cartesian grid. A body-fitted solver puts grids around the body instead: curvilinear structured blocks, overlapping each other and a background. This chapter builds an O-grid around an elliptic cylinder and checks it, locates points in it and interpolates from it (steps 1 and 2), then assembles it with a Cartesian background: holes, donors, receptors, fringes and the connectivity file a solver reads (step 3), the orphans of a too-narrow overlap (step 4), and the assembly of a moving O-grid (step 5). Everything is 2-D; the 3-D types and routines have the same names without _2d.

An O-grid

use fossil,     only : curvilinear_block_2d_t, topo_record_t, detect_topology_2d, TOPO_ABUT, TOPO_BOUNDARY, &
                       CELL_VALID, CELL_INVALID
use penf,       only : I4P, R8P

implicit none
integer(I4P), parameter      :: NI = 49, NJ = 17            ! O-grid vertices: around the body, out from the wall
real(R8P),    parameter      :: PI = acos(-1._R8P)
real(R8P),    parameter      :: A = 1._R8P, B = 0.4_R8P      ! semi-axes of the elliptic cylinder
type(curvilinear_block_2d_t) :: ogrid(1), folded       ! a set of one block
type(topo_record_t)          :: rec
real(R8P)                    :: u(NI, NJ), v(NI, NJ), c, mu0, mu, s
integer(I4P)                 :: i, j, r

c = sqrt(A**2 - B**2) ; mu0 = atanh(B / A)                   ! elliptic coordinates: the wall is mu = mu0
do j = 1, NJ
   s = real(j - 1, R8P) / (NJ - 1)
   mu = mu0 + 0.78_R8P * (exp(2._R8P * s) - 1._R8P) / (exp(2._R8P) - 1._R8P)   ! layers clustered at the wall
   do i = 1, NI
      u(i, j) = c * cosh(mu) * cos(2 * PI * (i - 1) / (NI - 1))
      v(i, j) = c * sinh(mu) * sin(2 * PI * (i - 1) / (NI - 1))
   enddo
enddo
u(NI,:) = u(1,:) ; v(NI,:) = v(1,:)                          ! the seam i = 1 = NI, written exactly
call ogrid(1)%initialize(u, v, name='ellipse')

• The O-grid has 49 × 17 vertices on the elliptic coordinates of an ellipse of semi-axes 1 and 0.4: i runs around the body, j out of it, with the layers clustered at the wall. The vertices are plain arrays u(ni, nj), v(ni, nj) (3-D: x, y, z(ni, nj, nk)); cell (i, j) is spanned by vertices i..i+1, j..j+1.
• initialize(u, v, name=) copies them. The last column repeats the first exactly: the seam.

A block alone knows nothing about its neighbours or its boundaries. Topology records say what each face is:

call detect_topology_2d(ogrid)                            ! the seam: one TOPO_ABUT record per side
do r = 1, ogrid(1)%nrec
   print '(A,I0,A,L1,A,I0,A,2I3)', 'record ', r, ': abut ', ogrid(1)%rec(r)%kind == TOPO_ABUT, ', face ', &
         ogrid(1)%rec(r)%face, ', transform', ogrid(1)%rec(r)%transform(1:2)
enddo
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [NI, 1, 1] ; rec%name = 'wall'
call ogrid(1)%add_record(rec)                                ! face 3 (-j): the wall
rec%face = 4 ; rec%lo = [1, NJ, 1] ; rec%hi = [NI, NJ, 1] ; rec%name = 'overlap'
call ogrid(1)%add_record(rec)                                ! face 4 (+j): the boundary that receives
print '(A,I0,A)', 'records after add_record: ', ogrid(1)%nrec, ' (2 abut, wall, overlap)'

• detect_topology_2d(blocks) proposes the records of a set of blocks from coincident vertices (exact by default, tolerance= for grids written with round-off). Here it finds the seam: one TOPO_ABUT record on face 1 (−i) and one on face 2 (+i), each pointing at the other, with the transform (1, 2): direction i of this side runs along direction i of the other. Walks and stencils cross the seam through these records. Periodic records come from translations=, polar axes as TOPO_SINGULAR.
• add_record appends records by hand. The physical boundaries are TOPO_BOUNDARY records with a name: the wall on face 3 (−j) and, on face 4 (+j), the boundary that receives from the background in the assembly. The assembly reads the names wall and overlap.

Checking the cells

A solver needs every cell's geometry and a guarantee that no cell is inverted:

call ogrid(1)%geometry(area=area)                            ! exact cell areas
print '(A,F10.6,A,F10.6)', 'O-grid area ', sum(area), ', annulus between the ellipses ', &
      PI * c**2 * (cosh(mu) * sinh(mu) - cosh(mu0) * sinh(mu0))
call ogrid(1)%validity(valid)                                ! exact: CELL_VALID, CELL_INVALID, CELL_SINGULAR
print '(A,I0,A,I0)', 'valid cells ', count(valid == CELL_VALID), ' of ', size(valid)
call ogrid(1)%quality(skewness=skew, orthogonality=ortho, aspect=aspect)
print '(A,F6.3,A,F6.3,A,F6.2,A,F5.2)', 'max skewness ', maxval(skew), ', min orthogonality ', minval(ortho), &
      ', aspect from ', maxval(aspect), ' at the wall to ', minval(aspect)
folded = ogrid(1)
folded%u(13, 5) = folded%u(13, 7) ; folded%v(13, 5) = folded%v(13, 7)  ! a vertex dragged two layers out
call folded%validity(valid)
print '(A,I0,A)', 'folded copy: ', count(valid == CELL_INVALID), ' invalid cells'
do j = 1, NJ - 1
   do i = 1, NI - 1
      if (valid(i, j) == CELL_INVALID) print '(A,2I3)', '  invalid cell', i, j
   enddo
enddo

• geometry(area=, centroid=, face_i=, face_j=) gives the exact cell geometry.
• validity(status) decides each cell exactly: CELL_VALID, CELL_INVALID, CELL_SINGULAR (a collapsed edge, a pole). In 2-D the edges are straight, so the four corner orientations (orient2d) decide; in 3-D the trilinear Jacobian is certified by its Bernstein coefficients.
• quality(skewness=, orthogonality=, aspect=, min_jacobian=) gives per-cell metrics: equiangle skewness, orthogonality of the faces to the line between neighbouring centroids, aspect ratio, and the certified minimum Jacobian over the mean.

$ tut13_a_block
record 1: abut T, face 1, transform  1  2
record 2: abut T, face 2, transform  1  2
records after add_record: 4 (2 abut, wall, overlap)
O-grid area  -5.992437, annulus between the ellipses   6.009585
valid cells 768 of 768
max skewness  0.103, min orthogonality  0.986, aspect from   8.06 at the wall to  1.23
folded copy: 2 invalid cells
  invalid cell 12  5
  invalid cell 13  5

• The two seam records come first; add_record brings the block to 4.
• The cells sum to −5.992: the sign is the block's handedness (i runs counter-clockwise, j outward, so the block is left-handed), not a defect. Its size is that of the polygon, 0.3 % below the area between the two ellipses, 6.010.
• All 768 cells are valid. The skewness is at most 0.103 and the orthogonality at least 0.986, since elliptic coordinates are orthogonal; the aspect ratio goes from 8.06 in the flat cells at the wall to 1.23 outside.
• The folded copy drags vertex (13, 5) two layers out, past the vertices of the layer above it: the two cells above it, (12, 5) and (13, 5), turn inside out and validity flags them; the two below are stretched but stay valid.

Locating points

Interpolation from the block starts by finding the cell that holds each point:

use fossil,     only : curvilinear_block_2d_t, detect_topology_2d, LOC_INSIDE, LOC_OUTSIDE, STENCIL_CENTRED
use penf,       only : I4P, R8P

implicit none
integer(I4P), parameter      :: NI = 49, NJ = 17, N = 200
real(R8P),    parameter      :: PI = acos(-1._R8P)
type(curvilinear_block_2d_t) :: ogrid(1)
real(R8P)                    :: u(NI, NJ), v(NI, NJ), pu(N), pv(N), s
real(R8P),    allocatable    :: xi(:,:), xi_walk(:,:)
integer(I4P), allocatable    :: cell(:,:), cell_walk(:,:), st(:), st_walk(:)
integer(I4P)                 :: i, fallbacks, jumps

call make_ogrid(u, v)                                        ! the O-grid of step 1
call ogrid(1)%initialize(u, v, name='ellipse')
call detect_topology_2d(ogrid)                               ! the seam records
call ogrid(1)%build_search                                   ! validity, handedness, the cell-box tree
do i = 1, N                                                  ! points on a spiral around the body
   s = real(i - 1, R8P) / (N - 1)
   pu(i) = (0.9_R8P + 0.9_R8P * s) * cos(14 * PI * s) ; pv(i) = (0.3_R8P + 1.1_R8P * s) * sin(14 * PI * s)
enddo
call ogrid(1)%locate_many(pu, pv, cell, xi, st)              ! robust: the tree, exact containment
print '(A,I0,A,I0,A,I0)', 'located ', count(st == LOC_INSIDE), ', outside ', count(st == LOC_OUTSIDE), ' of ', N
print '(A,I0,A,2I3,A,2F7.4)', 'point ', 100, ': cell', cell(:, 100), ', xi', xi(:, 100)
call ogrid(1)%locate_many(pu, pv, cell_walk, xi_walk, st_walk, walk=.true., jumps=jumps, fallbacks=fallbacks)
print '(A,L1,A,I0,A,I0)', 'walk equal to the robust path: ', all(cell_walk == cell) .and. all(xi_walk == xi) .and. &
      all(st_walk == st), ', cells crossed ', jumps, ', fallbacks ', fallbacks

• build_search stores the cell validity and the handedness and builds a bounding-box tree over the cells.
• locate_many(pu, pv, cell, xi, status) returns per point the cell (i, j), the local coordinates ξ in [0, 1]² and a status: LOC_INSIDE, LOC_OUTSIDE, or, when only a cell that is not valid contains the point, LOC_WARPED or LOC_DEGENERATE. Containment is exact in 2-D; a point on a face shared by two cells belongs to the low one.
• walk=.true. jumps from the previous point's cell towards the point (stencil jumping): cheap when consecutive points are close, as the vertices of a grid line are. It hands the point to the tree when it cannot settle it, for instance when the walk reaches the wall or the outer boundary; jumps= counts the cells crossed, fallbacks= the points handed over. The result equals the robust one bit for bit.

Interpolating with stencils

real(R8P),    allocatable    :: w(:,:)
integer(I4P), allocatable    :: nodes(:,:,:), sst(:)
real(R8P)                    :: err
integer(I4P)                 :: q, k
do q = 2, 4
   call ogrid(1)%stencils(cell, xi, q, nodes, w, sst)        ! q x q Lagrange stencils on the vertices
   err = 0._R8P
   do i = 1, N
      if (st(i) /= LOC_INSIDE) cycle
      err = max(err, abs(sum([(w(k, i) * fun(u(nodes(1, k, i), nodes(2, k, i)), v(nodes(1, k, i), nodes(2, k, i))), &
                                k = 1, q**2)]) - fun(pu(i), pv(i))))
   enddo
   print '(A,I0,A,I0,A,ES9.2)', 'q = ', q, ': centred stencils ', count(sst == STENCIL_CENTRED .and. st == LOC_INSIDE), &
         ', max error of sin(2x) cos(y) ', err
enddo

• stencils(cell, xi, q, nodes, weights, status) builds a q × q Lagrange stencil on the vertices around each located point: nodes(:, k, n) are vertex indices, weights(k, n) their weights. Near the wall and the outer boundary it is shifted one-sided (STENCIL_SHIFTED); across the seam it takes the vertices of the other side through the records.
• For q > 2, ξ is recomputed on the isoparametric map of degree q − 1 of the stencil's own vertices (order_matched=, on by default): the bilinear ξ alone would cap the interpolation at second order.

$ tut13_b_locate
located 165, outside 35 of 200
point 100: cell 24 14, xi 0.2754 0.1635
walk equal to the robust path: T, cells crossed 314, fallbacks 38
q = 2: centred stencils 165, max error of sin(2x) cos(y)  9.78E-03
q = 3: centred stencils 149, max error of sin(2x) cos(y)  1.75E-03
q = 4: centred stencils 133, max error of sin(2x) cos(y)  8.53E-04

• 165 of the 200 points are inside the block. Point 100 lies in cell (24, 14) at ξ = (0.275, 0.164).
• The walk returns the same cells and ξ bit for bit; it crosses 314 cells and hands 38 points to the tree, those whose walk ends at the wall or the outer boundary.
• The interpolation error of sin(2x) cos(y) falls from 9.78e-3 at q = 2 to 1.75e-3 at q = 3 and 8.53e-4 at q = 4. Fewer stencils are centred as q grows (165, 149, 133): the others are shifted at the wall or the outer boundary.

Cutting holes

The O-grid now overlaps a Cartesian background of 41 × 33 vertices at spacing 0.15. The background covers the body too: its vertices inside the body must be blanked.

use fossil,     only : curvilinear_block_2d_t, topo_record_t, detect_topology_2d, TOPO_BOUNDARY
use fossil,     only : curve_object, wall_cutter_2d, hole_flags_2d, hole_map_2d, block_logical_t, donor_set_t, &
                       donor_candidates_2d, overset_mask_t, classify_2d, fringe_levels_2d, overlap_width_2d, &
                       connectivity_t, build_connectivity, save_connectivity, SELECT_EXPLICIT, FRINGE_INNER_1, &
                       FRINGE_INNER_2, FRINGE_OUTER_1, FRINGE_OUTER_2
use penf,       only : I4P, R8P

implicit none
integer(I4P), parameter            :: NI = 49, NJ = 17, MI = 41, MJ = 33
type(curvilinear_block_2d_t)       :: blocks(2)              ! 1: the O-grid, 2: the background
type(topo_record_t)                :: rec
type(curve_object)                 :: cutter
type(block_logical_t), allocatable :: holes(:), holes_map(:)
type(donor_set_t),     allocatable :: donors(:)
type(overset_mask_t),  allocatable :: masks(:)
type(connectivity_t)               :: conn
real(R8P)                          :: u(NI, NJ), v(NI, NJ), x(MI, MJ), y(MI, MJ)
real(R8P),             allocatable :: d(:,:)
integer(I4P)                       :: i, j, b, iterations, circular, un
character(200)                     :: line

call make_ogrid(u, v)                                        ! the O-grid of step 1, with its records
call blocks(1)%initialize(u, v, name='ellipse')
call detect_topology_2d(blocks(1:1))
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [NI, 1, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec)
rec%face = 4 ; rec%lo = [1, NJ, 1] ; rec%hi = [NI, NJ, 1] ; rec%name = 'overlap'
call blocks(1)%add_record(rec)
do j = 1, MJ                                                 ! a Cartesian background, spacing 0.15
   do i = 1, MI
      x(i, j) = -3._R8P + 0.15_R8P * (i - 1) ; y(i, j) = -2.4_R8P + 0.15_R8P * (j - 1)
   enddo
enddo
call blocks(2)%initialize(x, y, name='background')
do b = 1, 2
   call blocks(b)%build_search
enddo

call wall_cutter_2d(blocks(1:1), cutter)                     ! the cutter: the edges of the faces named 'wall'
print '(A,I0,A)', 'cutter: ', cutter%get_vertices_number(), ' vertices, one closed loop'
call hole_flags_2d(blocks, cutter, holes, cuts=[.false., .true.]) ! the background only: the O-grid's wall is the cutter
print '(A,I0,A,I0)', 'background holes ', count(holes(2)%v), ' of ', size(holes(2)%v)
call hole_map_2d(blocks, cutter, holes_map, [64, 64], cuts=[.false., .true.])  ! a Cartesian hole map instead
print '(A,L1)', 'hole map gives the same holes: ', all(holes_map(2)%v .eqv. holes(2)%v)

• wall_cutter_2d(blocks, cutter) builds the cutter, a curve_object, from the edges of the faces recorded as wall, oriented with the solid on the left. Several walls make several loops; overlapping wall patches are trimmed.
• hole_flags_2d(blocks, cutter, holes, cuts=) flags the vertices inside the cutter by its fast winding number: exact on a closed wall. cuts= says which blocks the cutter cuts: the background, not the O-grid, whose wall lies on the cutter. offset= blanks below that signed distance instead of 0, a band off the wall.
• hole_map_2d(blocks, cutter, holes, n) is the alternative of the founding codes: a Cartesian map of n(1) × n(2) elements over the cutter, the elements the wall meets found exactly and the others painted outside from the corners. A vertex in an uncut element takes its status; one in a cut element looks along a grid axis to the first uncut element and counts the exact wall crossings in between. Here it blanks the same vertices. xray_flags_2d is the X-ray variant, exact only on its rays.

Donors and receptors

The vertices next to the holes and on the O-grid's outer boundary have no equation of their own: they receive values interpolated from the other grid.

call donor_candidates_2d(blocks, donors)                     ! every vertex located in the other block
print '(A,I0,A,I0)', 'donor candidates: O-grid ', size(donors(1)%block), ', background ', size(donors(2)%block)
call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes, iterations=iterations, &
                 circular=circular)                          ! 3 x 3 donor stencils, explicit interpolation
do b = 1, 2
   print '(2X,A,A,I0,A,I0,A,I0,A,I0)', blocks(b)%name(1:10), ': field ', count(masks(b)%iota == 1), ', holes ', &
         count(masks(b)%iota == 0), ', receptors ', count(masks(b)%iota < 0), ', orphans ', count(masks(b)%orphan > 0)
enddo
print '(A,I0,A,I0)', 'fixed point in ', iterations, ' iterations, circular pairs ', circular

• donor_candidates_2d(blocks, donors) locates every vertex of every block in every other block and keeps all candidates, each with its donor cell, ξ, the donor and receptor cell sizes and the certified minimum Jacobian of the donor.
• classify_2d(blocks, donors, masks, SELECT_EXPLICIT, q, holes=) gives the mask per vertex: iota = 1 for a field point, 0 for a hole, −k for a receptor from block k; donor is its chosen candidate. With the explicit rule the ring around the holes and the vertices of the overlap faces must receive; a candidate is valid when its q × q donor stencil holds no hole and no receptor, and candidates are tried by increasing donor size. Validity depends on the neighbours, so the classification is iterated to a fixed point; circular= counts receptor pairs in each other's stencils. SELECT_IMPLICIT chooses by donor cell size instead: a vertex receives when a smaller donor cell exists.

Fringes and connectivity

call fringe_levels_2d(blocks, masks)                         ! INNER 1, 2 around the holes; OUTER 1, 2 at the overlap
print '(A,4(1X,I0))', 'fringe INNER 1, INNER 2, OUTER 1, OUTER 2:', &
      sum([(count(masks(b)%fringe == FRINGE_INNER_1), b = 1, 2)]), sum([(count(masks(b)%fringe == FRINGE_INNER_2), b = 1, 2)]), &
      sum([(count(masks(b)%fringe == FRINGE_OUTER_1), b = 1, 2)]), sum([(count(masks(b)%fringe == FRINGE_OUTER_2), b = 1, 2)])
call overlap_width_2d(blocks, masks, d)                      ! Chesshire's width between the two fringes
print '(A,F7.4)', 'overlap width ', d(1, 2)
call build_connectivity(masks, donors, 3_I4P, 2_I4P, conn)   ! the receptor list a solver reads
call save_connectivity(conn, 'tut13-connectivity.txt')
open(newunit=un, file='tut13-connectivity.txt', action='read')
do i = 1, 5                                                  ! the header and the first receptor
   read(un, '(A)') line
   print '(A)', trim(line)
enddo
close(un)

• fringe_levels_2d sets masks(b)%fringe: INNER 1 and 2 around the holes, OUTER 1 and 2 at the overlap boundary, the levels a solver with a wider stencil needs.
• overlap_width_2d is Chesshire's width, the least distance between the receptors of the two blocks.
• build_connectivity(masks, donors, q, 2, conn) makes the receptor list and save_connectivity(conn, file) writes it as text: per block the counts by class, then one line per receptor with its block and vertex, the donor block, the donor cell, ξ, the stencil width, a low-order flag and the donor-to-receptor size ratio.

$ tut13_c_assembly
cutter: 48 vertices, one closed loop
background holes 57 of 1353
hole map gives the same holes: T
donor candidates: O-grid 833, background 268
  ellipse: field 784, holes 0, receptors 49, orphans 0
  background: field 1256, holes 57, receptors 40, orphans 0
fixed point in 2 iterations, circular pairs 0
fringe INNER 1, INNER 2, OUTER 1, OUTER 2: 40 48 49 49
overlap width  0.5401
# FOSSIL overset connectivity (issue #72)
# block 1 field 784 holes 0 receptors 49 orphans 0
# block 2 field 1256 holes 57 receptors 40 orphans 0
# receptors 89
1 785 2 32 16 1   9.7015928268724919E-02   1.0000000000000000E+00   0.0000000000000000E+00 3 3 1 F   9.5375240008439077E-01

• The cutter has 48 vertices, one per wall edge; it blanks 57 background vertices, and the hole map agrees.
• The 40 background vertices around the holes receive from the O-grid and the 49 vertices of the O-grid's outer boundary from the background. No orphan, two iterations of the fixed point, no circular pair; the overlap width is 0.540.
• The connectivity file lists the 89 receptors; the first is the outer vertex 785 of the O-grid, (1, 17), in background cell (32, 16) at ξ = (0.097, 1.0): on that cell's top edge, which the low-side rule gives to the lower cell.

Orphans

The overlap has to be wide enough for the stencils. The same assembly with the O-grid cut to fewer layers:

use fossil,     only : curvilinear_block_2d_t, topo_record_t, detect_topology_2d, TOPO_BOUNDARY, curve_object, &
                       wall_cutter_2d, hole_flags_2d, block_logical_t, donor_set_t, donor_candidates_2d, &
                       overset_mask_t, classify_2d, SELECT_EXPLICIT, ORPHAN_STENCIL_HOLE, ORPHAN_RECEPTOR_DONOR
use penf,       only : I4P, R8P

implicit none
integer(I4P), parameter            :: NI = 49, NJ = 17, MI = 41, MJ = 33
type(curvilinear_block_2d_t)       :: blocks(2)
type(topo_record_t)                :: rec
type(curve_object)                 :: cutter
type(block_logical_t), allocatable :: holes(:)
type(donor_set_t),     allocatable :: donors(:)
type(overset_mask_t),  allocatable :: masks(:)
real(R8P)                          :: u(NI, NJ), v(NI, NJ), x(MI, MJ), y(MI, MJ)
integer(I4P)                       :: i, j, k

call make_ogrid(u, v)                                        ! the O-grid of step 1
do j = 1, MJ                                                 ! the background of step 3
   do i = 1, MI
      x(i, j) = -3._R8P + 0.15_R8P * (i - 1) ; y(i, j) = -2.4_R8P + 0.15_R8P * (j - 1)
   enddo
enddo
call blocks(2)%initialize(x, y, name='background')
call blocks(2)%build_search
do k = 10, 16, 2                                             ! the first k - 1 cell layers of the O-grid only
   call blocks(1)%destroy
   call blocks(1)%initialize(u(:, 1:k), v(:, 1:k), name='thin')
   call detect_topology_2d(blocks(1:1))
   rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [NI, 1, 1] ; rec%name = 'wall'
   call blocks(1)%add_record(rec)
   rec%face = 4 ; rec%lo = [1, k, 1] ; rec%hi = [NI, k, 1] ; rec%name = 'overlap'
   call blocks(1)%add_record(rec)
   call blocks(1)%build_search
   call wall_cutter_2d(blocks(1:1), cutter)
   call hole_flags_2d(blocks, cutter, holes, cuts=[.false., .true.])
   call donor_candidates_2d(blocks, donors)
   call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes)
   print '(A,I2,A,I2,A,I2,A,I2,A,I2,A)', 'O-grid of ', k - 1, ' layers: receptors ', count(masks(1)%iota < 0), &
         ', orphans ', count(masks(1)%orphan > 0), ' (stencil on a hole ', &
         count(masks(1)%orphan == ORPHAN_STENCIL_HOLE), ', receptor in the stencil ', &
         count(masks(1)%orphan == ORPHAN_RECEPTOR_DONOR), ')'
enddo

• A receptor without a valid donor is an orphan; masks(b)%orphan gives the cause: ORPHAN_STENCIL_HOLE (every candidate stencil touches a hole), ORPHAN_RECEPTOR_DONOR (it holds a receptor, explicit interpolation), ORPHAN_NO_CANDIDATE (no donor cell at all), ORPHAN_BOUNDARY, ORPHAN_WALL.

$ tut13_d_orphans
O-grid of  9 layers: receptors  0, orphans 49 (stencil on a hole 16, receptor in the stencil 33)
O-grid of 11 layers: receptors 18, orphans 31 (stencil on a hole  4, receptor in the stencil 27)
O-grid of 13 layers: receptors 41, orphans  8 (stencil on a hole  0, receptor in the stencil  8)
O-grid of 15 layers: receptors 49, orphans  0 (stencil on a hole  0, receptor in the stencil  0)

• With 15 layers all 49 outer vertices find a donor; with 13 layers 8 are orphans, with 11 layers 31, with 9 all of them. The O-grid is thinnest at the ends of the long axis, and the orphans start there: their stencils reach the receptors around the holes, or the holes themselves. A solver cannot interpolate there; the grid has to change.

A moving near-body grid

When the body moves, its O-grid moves with it and the assembly is redone at every step. The O-grid does not need to be rewritten: it stays in its body frame, with its search tree and its cutter built once, and the kernels take its pose.

use fossil,     only : curvilinear_block_2d_t, topo_record_t, detect_topology_2d, TOPO_BOUNDARY, rigid_pose_2d_t, &
                       pose_2d, to_world, curve_object, wall_cutter_2d, hole_flags_2d, block_logical_t, donor_set_t, &
                       donor_candidates_2d, overset_mask_t, classify_2d, exposed_points, SELECT_EXPLICIT
use penf,       only : I4P, R8P

implicit none
integer(I4P), parameter            :: NI = 49, NJ = 17, MI = 41, MJ = 33, NF = 12   ! NF steps per period
real(R8P),    parameter            :: PI = acos(-1._R8P)
type(curvilinear_block_2d_t)       :: blocks(2)
type(topo_record_t)                :: rec
type(curve_object)                 :: cutter
type(rigid_pose_2d_t)              :: poses(2)
type(block_logical_t), allocatable :: holes(:), exposed(:)
type(donor_set_t),     allocatable :: donors(:)
type(overset_mask_t),  allocatable :: masks(:), old(:)
real(R8P)                          :: u(NI, NJ), v(NI, NJ), x(MI, MJ), y(MI, MJ), angle
integer(I4P)                       :: i, j, k, violations, n

call make_ogrid(u, v)                                        ! the O-grid of step 1, in its body frame
call blocks(1)%initialize(u, v, name='ellipse')
call detect_topology_2d(blocks(1:1))
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [NI, 1, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec)
rec%face = 4 ; rec%lo = [1, NJ, 1] ; rec%hi = [NI, NJ, 1] ; rec%name = 'overlap'
call blocks(1)%add_record(rec)
do j = 1, MJ
   do i = 1, MI
      x(i, j) = -3._R8P + 0.15_R8P * (i - 1) ; y(i, j) = -2.4_R8P + 0.15_R8P * (j - 1)
   enddo
enddo
call blocks(2)%initialize(x, y, name='background')
call blocks(1)%build_search ; call blocks(2)%build_search    ! once: the O-grid moves by its pose only
call wall_cutter_2d(blocks(1:1), cutter)                     ! once, in the body frame
print '(A)', 'step  pitch  holes  receptors (bg + O-grid)  orphans  exposed  holes to field'
do k = 0, NF
   angle = PI / 6._R8P * sin(2 * PI * mod(k, NF) / NF)     ! pitch 30 degrees, heave 0.3
   poses(1) = pose_2d(angle=angle, translation=[0._R8P, 0.3_R8P * sin(2 * PI * k / NF + PI / 2._R8P)])
   poses(2) = pose_2d(angle=0._R8P)                          ! the background stays
   call hole_flags_2d(blocks, cutter, holes, cuts=[.false., .true.], poses=poses, cutter_pose=poses(1))
   call donor_candidates_2d(blocks, donors, poses=poses)
   call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes)
   if (k == 0) old = masks
   call exposed_points(old, masks, exposed, violations)      ! uncovered since the previous step
   print '(I4,F7.1,I7,I8,A,I3,I15,I9,I12)', k, angle * 180._R8P / PI, count(masks(2)%iota == 0), &
         count(masks(2)%iota < 0), ' +', count(masks(1)%iota < 0), count(masks(1)%orphan > 0) + count(masks(2)%orphan > 0), &
         count(exposed(2)%v), violations
   old = masks
enddo

• poses= on hole_flags_2d and donor_candidates_2d gives each block's pose (a rigid_pose_2d_t, chapter 11); cutter_pose= the cutter's. The background takes the identity.
• exposed_points(old, new, exposed, violations) compares two classifications: exposed flags the vertices that were holes or receptors and are now field points; they need values before the solver advances them. violations counts the holes that became field points directly, without receiving first (Meakin & Suhs): the step was too long for the fringe.

n = 0
do k = 0, 2 * NF                                             ! the same period in steps half as long
   poses(1) = pose_2d(angle=PI / 6._R8P * sin(PI * mod(k, 2 * NF) / NF), &
                      translation=[0._R8P, 0.3_R8P * sin(PI * k / NF + PI / 2._R8P)])
   call hole_flags_2d(blocks, cutter, holes, cuts=[.false., .true.], poses=poses, cutter_pose=poses(1))
   call donor_candidates_2d(blocks, donors, poses=poses)
   call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes)
   if (k > 0) call exposed_points(old, masks, exposed, violations)
   if (k > 0) n = n + violations
   old = masks
enddo
print '(A,I0,A,I0)', 'steps per period ', 2 * NF, ': holes to field ', n

$ tut13_e_moving
step  pitch  holes  receptors (bg + O-grid)  orphans  exposed  holes to field
   0    0.0     57      40 + 49              0        0           0
   1   15.0     56      42 + 49              0       11           0
   2   26.0     53      44 + 49              0       10           0
   3   30.0     53      44 + 49              0       16           0
   4   26.0     53      44 + 49              0       16           1
   5   15.0     56      42 + 49              0        9           0
   6    0.0     57      40 + 49              0       12           0
   7  -15.0     56      42 + 49              0       11           0
   8  -26.0     53      44 + 49              0       10           0
   9  -30.0     53      44 + 49              0       16           0
  10  -26.0     53      44 + 49              0       16           1
  11  -15.0     56      42 + 49              0        9           0
  12    0.0     57      40 + 49              0       12           0
steps per period 24: holes to field 0

• Over one period of 12 steps (pitch ±30 degrees, heave ±0.3) the holes go from 57 to 53 and back, the background receptors from 40 to 44; no orphan at any pose. After the first step, 9 to 16 background vertices are exposed per step.
• At steps 4 and 10, the first steps back from the extreme pitch, one hole becomes a field point directly: the step is too long for the fringe there. With 24 steps per period no hole does.

What you learned

A curvilinear_block_2d_t with its topology records (detect_topology_2d, add_record), exact cell geometry, validity and quality, located points (build_search, locate_many, walk=) and Lagrange stencils (stencils). An overset assembly from kernels: wall_cutter_2d, hole_flags_2d (or hole_map_2d), donor_candidates_2d, classify_2d, fringe_levels_2d, orphans by cause, build_connectivity and save_connectivity; a moving block by poses=, and exposed_points between two steps. Reference: Curvilinear blocks, Overset assembly.

Next: 14. Contact: distance and collision.