/**
 * embedded3 R6-1: the body's φ and the cut geometry on the device — the
 * host `cut.rs::CutGeometry::build_from` expression for expression in fp64
 * (compiled with FMA contraction off; `/` and `sqrt` are IEEE round-to-
 * nearest in double on the device, as on the host), and the flag test's
 * φ (`tests/embedded3_flag_wake.rs`: the circle, the capsule around the
 * centreline polyline, the span cut with rounded edges, the fillet union).
 *
 * Kernels:
 *   e3_geom_phi    per corner: the narrow band's keep-or-evaluate (bound −
 *                  motion > band keeps φ), φ from the polyline otherwise;
 *   e3_geom_faces  per face of one component: the aperture by the two
 *                  Kuhn triangles and the face-centre φ, where a corner was
 *                  evaluated (every face on a full pass);
 *   e3_geom_cells  per cell: the fluid volume by the six Kuhn tetrahedra
 *                  and the wall vector by closure of the face apertures;
 *   e3_geom_gather compact copies for the host mirror.
 */

struct GeomSdf {
    double cx, cy, rc;      /* the circle */
    double zc, span, r_edge;/* the span cut */
    double half;            /* capsule half-thickness */
    double fillet;          /* root fillet radius (0: min) */
    int cyl_cut, flag_cut;  /* cut to the span */
    int npts;               /* polyline points (x, y interleaved) */
    int nst, ns;            /* R8-c plate: stations and points per station (nst 0: the polyline) */
    int flat_tip;           /* R8-h: 1 = the flat tip (corners rounded to tip_corner); 0 = the capsule */
    double tip_corner;
};

struct GeomGrid {
    int nx, ny, nz;
    double dx, dy, dz;
};

__device__ __forceinline__ double rs_min(double a, double b) { return b < a ? b : a; }
__device__ __forceinline__ double rs_max(double a, double b) { return b > a ? b : a; }

/* The span cut with rounded edges (flag_3d / cylinder_3d). */
__device__ __forceinline__ double span_cut(double d2, double z, const GeomSdf& s)
{
    double r = s.r_edge;
    double q1 = d2 + r;
    double q2 = fabs(z - s.zc) - 0.5 * s.span + r;
    double m1 = rs_max(q1, 0.0), m2 = rs_max(q2, 0.0);
    double outside = sqrt(m1 * m1 + m2 * m2);
    return outside + rs_min(rs_max(q1, q2), 0.0) - r;
}

/* R8-h: the tip's axial distance along the last segment a→b from b (plate.rs tip_axial). */
__device__ __forceinline__ double tip_axial(double x, double y, double ax, double ay, double bx, double by)
{
    double ex = bx - ax, ey = by - ay;
    double l = sqrt(ex * ex + ey * ey);
    return ((x - bx) * ex + (y - by) * ey) / l;
}

/* R8-h: the flat tip's rounded cap (plate.rs flat_cap). */
__device__ __forceinline__ double flat_cap(double dl, double a, double r)
{
    double q1 = dl + r;
    double q2 = a + r;
    double m1 = rs_max(q1, 0.0), m2 = rs_max(q2, 0.0);
    double outside = sqrt(m1 * m1 + m2 * m2);
    return outside + rs_min(rs_max(q1, q2), 0.0) - r;
}

/*
 * R8-c: the deformed plate (`plate.rs`, the host twin expression for
 * expression): `poly` holds the stations' rows (x, y interleaved, row-major
 * by station) followed by the stations' z. The stations bracketing z are
 * interpolated linearly into one polyline (outside their range: the end
 * station as it is), the in-plane closest point taken as the polyline's,
 * the distance corrected for the spanwise slope d / sqrt(1 + (e·c)²);
 * with `vel`, the velocity at the closest point the same way.
 */
__device__ double plate_dist(double x, double y, double z, const GeomSdf& s,
                             const double* __restrict__ P, const double* __restrict__ V,
                             double* vx, double* vy, double* axial)
{
    int nst = s.nst, ns = s.ns;
    const double* Z = P + 2 * (long long) nst * ns;
    int k = 0, interp = 0;
    double w = 0.0;
    if (nst > 1) {
        if (z <= Z[0]) {
            k = 0;
        } else if (z >= Z[nst - 1]) {
            k = nst - 2;
        } else {
            while (k + 1 < nst && Z[k + 1] <= z) ++k;
            if (k > nst - 2) k = nst - 2;
        }
        interp = 1;
        w = (z - Z[k]) / (Z[k + 1] - Z[k]);
        /* beyond the end stations: at most half an interval extrapolated */
        if (w < -0.5) w = -0.5;
        if (w > 1.5) w = 1.5;
    }
    const double* R0 = P + 2 * (long long) k * ns;
    const double* R1 = interp ? R0 + 2 * ns : R0;
    double best = 1.0 / 0.0, ub = 0.0, qbx = 0.0, qby = 0.0;
    int mb = 0;
    for (int m = 0; m + 1 < ns; ++m) {
        double ax, ay, bx, by;
        if (interp) {
            ax = R0[2 * m] + w * (R1[2 * m] - R0[2 * m]);
            ay = R0[2 * m + 1] + w * (R1[2 * m + 1] - R0[2 * m + 1]);
            bx = R0[2 * m + 2] + w * (R1[2 * m + 2] - R0[2 * m + 2]);
            by = R0[2 * m + 3] + w * (R1[2 * m + 3] - R0[2 * m + 3]);
        } else {
            ax = R0[2 * m]; ay = R0[2 * m + 1];
            bx = R0[2 * m + 2]; by = R0[2 * m + 3];
        }
        double ex = bx - ax, ey = by - ay;
        double l2 = ex * ex + ey * ey;
        double u = ((x - ax) * ex + (y - ay) * ey) / l2;
        if (u < 0.0) u = 0.0;
        if (u > 1.0 && !(s.flat_tip && m + 2 == ns)) u = 1.0;
        double px = ax + u * ex, py = ay + u * ey;
        double qx = x - px, qy = y - py;
        double d = sqrt(qx * qx + qy * qy);
        if (d < best) {
            best = d;
            mb = m;
            ub = u > 1.0 ? 1.0 : u;
            qbx = qx;
            qby = qy;
        }
    }
    if (interp && best > 0.0) {
        double dz = Z[k + 1] - Z[k];
        double cax = R1[2 * mb] - R0[2 * mb], cay = R1[2 * mb + 1] - R0[2 * mb + 1];
        double cbx = R1[2 * mb + 2] - R0[2 * mb + 2], cby = R1[2 * mb + 3] - R0[2 * mb + 3];
        double cx = (cax + ub * (cbx - cax)) / dz;
        double cy = (cay + ub * (cby - cay)) / dz;
        double qn = (qbx * cx + qby * cy) / best;
        best = best / sqrt(1.0 + qn * qn);
    }
    if (s.flat_tip && axial) {
        int m = ns - 2;
        double ax, ay, bx, by;
        if (interp) {
            ax = R0[2 * m] + w * (R1[2 * m] - R0[2 * m]);
            ay = R0[2 * m + 1] + w * (R1[2 * m + 1] - R0[2 * m + 1]);
            bx = R0[2 * m + 2] + w * (R1[2 * m + 2] - R0[2 * m + 2]);
            by = R0[2 * m + 3] + w * (R1[2 * m + 3] - R0[2 * m + 3]);
        } else {
            ax = R0[2 * m]; ay = R0[2 * m + 1];
            bx = R0[2 * m + 2]; by = R0[2 * m + 3];
        }
        *axial = tip_axial(x, y, ax, ay, bx, by);
    }
    if (V) {
        const double* V0 = V + 2 * (long long) k * ns;
        const double* V1 = interp ? V0 + 2 * ns : V0;
        double avx, avy, bvx, bvy;
        if (interp) {
            avx = V0[2 * mb] + w * (V1[2 * mb] - V0[2 * mb]);
            avy = V0[2 * mb + 1] + w * (V1[2 * mb + 1] - V0[2 * mb + 1]);
            bvx = V0[2 * mb + 2] + w * (V1[2 * mb + 2] - V0[2 * mb + 2]);
            bvy = V0[2 * mb + 3] + w * (V1[2 * mb + 3] - V0[2 * mb + 3]);
        } else {
            avx = V0[2 * mb]; avy = V0[2 * mb + 1];
            bvx = V0[2 * mb + 2]; bvy = V0[2 * mb + 3];
        }
        *vx = avx + ub * (bvx - avx);
        *vy = avy + ub * (bvy - avy);
    }
    return best;
}

__device__ double geom_phi_at(double x, double y, double z, const GeomSdf& s, const double* __restrict__ poly)
{
    /* the circle */
    double ex0 = x - s.cx, ey0 = y - s.cy;
    double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
    if (s.cyl_cut) dc = span_cut(dc, z, s);
    /* the capsule: distance to the polyline (or the plate, R8-c) */
    double best = 1.0 / 0.0, axial = 0.0;
    if (s.nst > 0) best = plate_dist(x, y, z, s, poly, nullptr, nullptr, nullptr, &axial);
    else {
        for (int m = 0; m + 1 < s.npts; ++m) {
            double ax = poly[2 * m], ay = poly[2 * m + 1];
            double bx = poly[2 * m + 2], by = poly[2 * m + 3];
            double ex = bx - ax, ey = by - ay;
            double l2 = ex * ex + ey * ey;
            double u = ((x - ax) * ex + (y - ay) * ey) / l2;
            if (u < 0.0) u = 0.0;
            if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0;
            double px = ax + u * ex, py = ay + u * ey;
            double qx = x - px, qy = y - py;
            double d = sqrt(qx * qx + qy * qy);
            if (d < best) best = d;
        }
        if (s.flat_tip) {
            int m = s.npts - 2;
            axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]);
        }
    }
    double df = best - s.half;
    if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner);
    if (s.flag_cut) df = span_cut(df, z, s);
    double r = s.fillet;
    if (r > 0.0 && dc < r && df < r) {
        double a = r - dc, b = r - df;
        return r - sqrt(a * a + b * b);
    }
    return rs_min(dc, df);
}

/* Per corner. `has_prev`: keep φ where the decayed bound stays beyond the band. */
extern "C" __global__ void e3_geom_phi(
    GeomGrid g, GeomSdf s, const double* __restrict__ poly,
    int has_prev, double band, double motion,
    double* __restrict__ phi, double* __restrict__ bound, unsigned char* __restrict__ touched)
{
    long long n = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    long long nn = (long long) (g.nx + 1) * (g.ny + 1) * (g.nz + 1);
    if (n >= nn) return;
    long long nxy = (long long) (g.nx + 1) * (g.ny + 1);
    int k = (int) (n / nxy);
    int j = (int) ((n / (g.nx + 1)) % (g.ny + 1));
    int i = (int) (n % (g.nx + 1));
    if (has_prev) {
        double b = bound[n] - motion;
        if (b > band) {
            bound[n] = b;
            touched[n] = 0;
            return;
        }
    }
    double v = geom_phi_at((double) i * g.dx, (double) j * g.dy, (double) k * g.dz, s, poly);
    phi[n] = v;
    bound[n] = fabs(v);
    touched[n] = 1;
}

__device__ __forceinline__ long long gnode(const GeomGrid& g, int k, int j, int i)
{
    return ((long long) k * (g.ny + 1) + j) * (g.nx + 1) + i;
}

/* cut.rs tri_area_fraction */
__device__ double tri_frac(double p0, double p1, double p2)
{
    double v[3] = {p0, p1, p2};
    int pos = (p0 >= 0.0) + (p1 >= 0.0) + (p2 >= 0.0);
    if (pos == 0) return 0.0;
    if (pos == 3) return 1.0;
    if (pos == 1) {
        int a = v[0] >= 0.0 ? 0 : (v[1] >= 0.0 ? 1 : 2);
        int b = (a + 1) % 3, c = (a + 2) % 3;
        return (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]));
    }
    int a = v[0] < 0.0 ? 0 : (v[1] < 0.0 ? 1 : 2);
    int b = (a + 1) % 3, c = (a + 2) % 3;
    return 1.0 - (v[a] / (v[a] - v[b])) * (v[a] / (v[a] - v[c]));
}

__device__ __forceinline__ double quad_frac(double q00, double q10, double q01, double q11)
{
    return 0.5 * (tri_frac(q00, q10, q11) + tri_frac(q00, q11, q01));
}

/* Per face of component c (0 u, 1 v, 2 w): aperture and face-centre φ. */
extern "C" __global__ void e3_geom_faces(
    GeomGrid g, int c, int full, const double* __restrict__ phi, const unsigned char* __restrict__ touched,
    double* __restrict__ a, double* __restrict__ d)
{
    long long f = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    int ni = g.nx + (c == 0), nj = g.ny + (c == 1), nk = g.nz + (c == 2);
    if (f >= (long long) ni * nj * nk) return;
    int k = (int) (f / ((long long) nj * ni));
    int j = (int) ((f / ni) % nj);
    int i = (int) (f % ni);
    long long n00, n10, n01, n11;
    if (c == 0) {        /* x-face: (j, k) */
        n00 = gnode(g, k, j, i); n10 = gnode(g, k, j + 1, i);
        n01 = gnode(g, k + 1, j, i); n11 = gnode(g, k + 1, j + 1, i);
    } else if (c == 1) { /* y-face: (i, k) */
        n00 = gnode(g, k, j, i); n10 = gnode(g, k, j, i + 1);
        n01 = gnode(g, k + 1, j, i); n11 = gnode(g, k + 1, j, i + 1);
    } else {             /* z-face: (i, j) */
        n00 = gnode(g, k, j, i); n10 = gnode(g, k, j, i + 1);
        n01 = gnode(g, k, j + 1, i); n11 = gnode(g, k, j + 1, i + 1);
    }
    if (!full && !(touched[n00] | touched[n10] | touched[n01] | touched[n11])) return;
    double q00 = phi[n00], q10 = phi[n10], q01 = phi[n01], q11 = phi[n11];
    a[f] = quad_frac(q00, q10, q01, q11);
    d[f] = 0.25 * (q00 + q10 + q01 + q11);
}

__device__ __forceinline__ double det3(const double* a, const double* b, const double* c)
{
    return a[0] * (b[1] * c[2] - b[2] * c[1]) - a[1] * (b[0] * c[2] - b[2] * c[0])
        + a[2] * (b[0] * c[1] - b[1] * c[0]);
}

__device__ double tet_vol(const double* p0, const double* p1, const double* p2, const double* p3)
{
    double e1[3] = {p1[0] - p0[0], p1[1] - p0[1], p1[2] - p0[2]};
    double e2[3] = {p2[0] - p0[0], p2[1] - p0[1], p2[2] - p0[2]};
    double e3[3] = {p3[0] - p0[0], p3[1] - p0[1], p3[2] - p0[2]};
    return fabs(det3(e1, e2, e3)) / 6.0;
}

__device__ __forceinline__ void lerp3(const double* a, const double* b, double t, double* out)
{
    out[0] = a[0] + t * (b[0] - a[0]);
    out[1] = a[1] + t * (b[1] - a[1]);
    out[2] = a[2] + t * (b[2] - a[2]);
}

/* cut.rs tet_fluid_volume */
__device__ double tet_fluid(const double pts[4][3], const double* v)
{
    double total = tet_vol(pts[0], pts[1], pts[2], pts[3]);
    int pos[4], neg[4], np = 0, nn = 0;
    for (int q = 0; q < 4; ++q) {
        if (v[q] >= 0.0) pos[np++] = q;
    }
    for (int q = 0; q < 4; ++q) {
        if (v[q] < 0.0) neg[nn++] = q;
    }
    if (np == 0) return 0.0;
    if (np == 4) return total;
    if (np == 1 || np == 3) {
        int a = np == 1 ? pos[0] : neg[0];
        const int* o = np == 1 ? neg : pos;
        double pb[3], pc[3], pd[3];
        lerp3(pts[a], pts[o[0]], v[a] / (v[a] - v[o[0]]), pb);
        lerp3(pts[a], pts[o[1]], v[a] / (v[a] - v[o[1]]), pc);
        lerp3(pts[a], pts[o[2]], v[a] / (v[a] - v[o[2]]), pd);
        double t = tet_vol(pts[a], pb, pc, pd);
        return np == 1 ? t : total - t;
    }
    int a = pos[0], b = pos[1], c = neg[0], d = neg[1];
    double pac[3], pad[3], pbc[3], pbd[3];
    lerp3(pts[a], pts[c], v[a] / (v[a] - v[c]), pac);
    lerp3(pts[a], pts[d], v[a] / (v[a] - v[d]), pad);
    lerp3(pts[b], pts[c], v[b] / (v[b] - v[c]), pbc);
    lerp3(pts[b], pts[d], v[b] / (v[b] - v[d]), pbd);
    return tet_vol(pts[a], pts[b], pbc, pbd) + tet_vol(pts[a], pac, pbc, pbd)
        + tet_vol(pts[a], pac, pad, pbd);
}

/* The Kuhn split around (0,0,0)–(1,1,1), corners as (x, y, z). */
__constant__ int KUHN[6][4][3] = {
    {{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {1, 1, 1}},
    {{0, 0, 0}, {1, 0, 0}, {1, 0, 1}, {1, 1, 1}},
    {{0, 0, 0}, {0, 1, 0}, {1, 1, 0}, {1, 1, 1}},
    {{0, 0, 0}, {0, 1, 0}, {0, 1, 1}, {1, 1, 1}},
    {{0, 0, 0}, {0, 0, 1}, {1, 0, 1}, {1, 1, 1}},
    {{0, 0, 0}, {0, 0, 1}, {0, 1, 1}, {1, 1, 1}},
};

/* Per cell: fluid volume fraction and the wall vector (3 interleaved). */
extern "C" __global__ void e3_geom_cells(
    GeomGrid g, int full, const double* __restrict__ phi, const unsigned char* __restrict__ touched,
    const double* __restrict__ a_u, const double* __restrict__ a_v, const double* __restrict__ a_w,
    double* __restrict__ vol, double* __restrict__ wall, unsigned char* __restrict__ touched_cell)
{
    long long idx = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    long long nxy = (long long) g.nx * g.ny;
    if (idx >= nxy * g.nz) return;
    int k = (int) (idx / nxy);
    int j = (int) ((idx / g.nx) % g.ny);
    int i = (int) (idx % g.nx);
    unsigned char t = 0;
    for (int dk = 0; dk < 2; ++dk)
        for (int dj = 0; dj < 2; ++dj)
            for (int di = 0; di < 2; ++di)
                t |= touched[gnode(g, k + dk, j + dj, i + di)];
    touched_cell[idx] = t;
    if (!full && !t) return;
    double fluid = 0.0;
    for (int q = 0; q < 6; ++q) {
        double pts[4][3];
        double vals[4];
        for (int m = 0; m < 4; ++m) {
            pts[m][0] = (double) KUHN[q][m][0];
            pts[m][1] = (double) KUHN[q][m][1];
            pts[m][2] = (double) KUHN[q][m][2];
            vals[m] = phi[gnode(g, k + KUHN[q][m][2], j + KUHN[q][m][1], i + KUHN[q][m][0])];
        }
        fluid += tet_fluid(pts, vals);
    }
    vol[idx] = fluid;
    double ax = g.dy * g.dz, ay = g.dx * g.dz, az = g.dx * g.dy;
    long long fu0 = ((long long) k * g.ny + j) * (g.nx + 1) + i;
    long long fv0 = ((long long) k * (g.ny + 1) + j) * g.nx + i;
    long long fw0 = ((long long) k * g.ny + j) * g.nx + i;
    double sx = (a_u[fu0 + 1] - a_u[fu0]) * ax;
    double sy = (a_v[fv0 + g.nx] - a_v[fv0]) * ay;
    double sz = (a_w[fw0 + nxy] - a_w[fw0]) * az;
    wall[3 * idx] = -sx;
    wall[3 * idx + 1] = -sy;
    wall[3 * idx + 2] = -sz;
}

/* out[w t + c] = src[w idx[t] + c], c < w (compact copies for the host mirror). */
extern "C" __global__ void e3_geom_gather(
    int n, int w, const unsigned int* __restrict__ idx, const double* __restrict__ src, double* __restrict__ out)
{
    int t = blockIdx.x * blockDim.x + threadIdx.x;
    if (t >= n) return;
    long long s = (long long) w * idx[t];
    for (int c = 0; c < w; ++c) out[(long long) w * t + c] = src[s + c];
}

/*
 * R6-2 step 2: the body's surface velocity on the device (the flag test's
 * host closure: the centreline's velocity interpolated at the capsule's
 * closest point where the span-cut capsule is not farther than the circle,
 * zero on the circle, no z component), and its two uses on a moving cut
 * mask:
 *   e3_geom_ub      the imposition band's faces: the velocity at the foot
 *                   (`Mask::surface_velocity_at`: the face centre plus its
 *                   centroid shift when the centroid foot is on, φ and its
 *                   unit gradient from the trilinear interpolant of the
 *                   corner φ, the foot x − s n) into the dense table and a
 *                   packed copy for the host mirror;
 *   e3_geom_impose  the band's solid faces (not open at the instant, off the
 *                   domain sides the host loops skip): the velocity at the
 *                   face centre into the field (`Mask::impose_from`, whose
 *                   ghost lists are empty on a cut mask);
 *   e3_geom_seam    the periodic seam: w at k = nz takes k = 0's where both
 *                   are solid.
 */
__device__ double body_velocity(double x, double y, double z, int c, const GeomSdf& s,
                                const double* __restrict__ poly, const double* __restrict__ vel)
{
    double best = 1.0 / 0.0, vx = 0.0, vy = 0.0, axial = 0.0;
    if (s.nst > 0) best = plate_dist(x, y, z, s, poly, vel, &vx, &vy, &axial);
    else {
        for (int m = 0; m + 1 < s.npts; ++m) {
            double ax = poly[2 * m], ay = poly[2 * m + 1];
            double bx = poly[2 * m + 2], by = poly[2 * m + 3];
            double ex = bx - ax, ey = by - ay;
            double l2 = ex * ex + ey * ey;
            double u = ((x - ax) * ex + (y - ay) * ey) / l2;
            if (u < 0.0) u = 0.0;
            if (u > 1.0 && !(s.flat_tip && m + 2 == s.npts)) u = 1.0;
            double px = ax + u * ex, py = ay + u * ey;
            double qx = x - px, qy = y - py;
            double d = sqrt(qx * qx + qy * qy);
            if (d < best) {
                best = d;
                double uv = u > 1.0 ? 1.0 : u;
                double avx = vel[2 * m], avy = vel[2 * m + 1];
                double bvx = vel[2 * m + 2], bvy = vel[2 * m + 3];
                vx = avx + uv * (bvx - avx);
                vy = avy + uv * (bvy - avy);
            }
        }
        if (s.flat_tip) {
            int m = s.npts - 2;
            axial = tip_axial(x, y, poly[2 * m], poly[2 * m + 1], poly[2 * m + 2], poly[2 * m + 3]);
        }
    }
    double df = best - s.half;
    if (s.flat_tip) df = flat_cap(df, axial, s.tip_corner);
    if (s.flag_cut) df = span_cut(df, z, s);
    double ex0 = x - s.cx, ey0 = y - s.cy;
    double dc = sqrt(ex0 * ex0 + ey0 * ey0) - s.rc;
    if (s.cyl_cut) dc = span_cut(dc, z, s);
    if (df <= dc) return c == 0 ? vx : (c == 1 ? vy : 0.0);
    return 0.0;
}

__device__ __forceinline__ double ub_node(const GeomGrid& g, const double* __restrict__ phi,
                                          long long i, long long j, long long k)
{
    if (i < 0) i = 0;
    if (i > g.nx) i = g.nx;
    if (j < 0) j = 0;
    if (j > g.ny) j = g.ny;
    if (k < 0) k = 0;
    if (k > g.nz) k = g.nz;
    return phi[(k * (g.ny + 1) + j) * (g.nx + 1) + i];
}

__device__ __forceinline__ double ub_lerp(double a, double b, double f) { return a + f * (b - a); }

/* The face centre of face `f` of component `c`. */
__device__ __forceinline__ void face_centre(const GeomGrid& g, int c, long long f, double* x)
{
    long long ni = g.nx + (c == 0), nj = g.ny + (c == 1);
    long long k = f / (nj * ni), j = (f / ni) % nj, i = f % ni;
    x[0] = ((double) i + (c == 0 ? 0.0 : 0.5)) * g.dx;
    x[1] = ((double) j + (c == 1 ? 0.0 : 0.5)) * g.dy;
    x[2] = ((double) k + (c == 2 ? 0.0 : 0.5)) * g.dz;
}

extern "C" __global__ void e3_geom_ub(
    GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel,
    const double* __restrict__ phi, int c, const unsigned int* __restrict__ faces, long long nf,
    const double* __restrict__ shift, int use_shift,
    double* __restrict__ ub, double* __restrict__ packed)
{
    long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    if (q >= nf) return;
    long long f = faces[q];
    double xc[3];
    face_centre(g, c, f, xc);
    double x0 = xc[0], x1 = xc[1], x2 = xc[2];
    if (use_shift) {
        x0 = x0 + shift[3 * f];
        x1 = x1 + shift[3 * f + 1];
        x2 = x2 + shift[3 * f + 2];
    }
    /* the interpolant's value and unit gradient */
    double gx = x0 / g.dx, gy = x1 / g.dy, gz = x2 / g.dz;
    long long i0 = (long long) floor(gx), j0 = (long long) floor(gy), k0 = (long long) floor(gz);
    double fx = gx - (double) i0, fy = gy - (double) j0, fz = gz - (double) k0;
    double c000 = ub_node(g, phi, i0, j0, k0), c100 = ub_node(g, phi, i0 + 1, j0, k0);
    double c010 = ub_node(g, phi, i0, j0 + 1, k0), c110 = ub_node(g, phi, i0 + 1, j0 + 1, k0);
    double c001 = ub_node(g, phi, i0, j0, k0 + 1), c101 = ub_node(g, phi, i0 + 1, j0, k0 + 1);
    double c011 = ub_node(g, phi, i0, j0 + 1, k0 + 1), c111 = ub_node(g, phi, i0 + 1, j0 + 1, k0 + 1);
    double c00 = ub_lerp(c000, c100, fx);
    double c10 = ub_lerp(c010, c110, fx);
    double c01 = ub_lerp(c001, c101, fx);
    double c11 = ub_lerp(c011, c111, fx);
    double c0 = ub_lerp(c00, c10, fy);
    double c1 = ub_lerp(c01, c11, fy);
    double sd = ub_lerp(c0, c1, fz);
    double dx0 = ub_lerp(c100 - c000, c110 - c010, fy);
    double dx1 = ub_lerp(c101 - c001, c111 - c011, fy);
    double px = ub_lerp(dx0, dx1, fz) / g.dx;
    double dy0 = ub_lerp(c010 - c000, c110 - c100, fx);
    double dy1 = ub_lerp(c011 - c001, c111 - c101, fx);
    double py = ub_lerp(dy0, dy1, fz) / g.dy;
    double dz0 = ub_lerp(c001 - c000, c101 - c100, fx);
    double dz1 = ub_lerp(c011 - c010, c111 - c110, fx);
    double pz = ub_lerp(dz0, dz1, fy) / g.dz;
    double norm = sqrt(px * px + py * py + pz * pz);
    double n0 = 1.0, n1 = 0.0, n2 = 0.0;
    if (norm > 0.0) {
        n0 = px / norm;
        n1 = py / norm;
        n2 = pz / norm;
    }
    double v = body_velocity(x0 - sd * n0, x1 - sd * n1, x2 - sd * n2, c, s, poly, vel);
    ub[f] = v;
    packed[q] = v;
}

extern "C" __global__ void e3_geom_impose(
    GeomGrid g, GeomSdf s, const double* __restrict__ poly, const double* __restrict__ vel,
    int c, const unsigned int* __restrict__ faces, long long nf,
    const int* __restrict__ open, double* __restrict__ field)
{
    long long q = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    if (q >= nf) return;
    long long f = faces[q];
    if (open[f] != 0) return;
    long long ni = g.nx + (c == 0), nj = g.ny + (c == 1);
    long long j = (f / ni) % nj, i = f % ni;
    /* the host loops: u over i in 1..nx, v over j in 1..ny, every w */
    if (c == 0 && (i == 0 || i == g.nx)) return;
    if (c == 1 && (j == 0 || j == g.ny)) return;
    double xc[3];
    face_centre(g, c, f, xc);
    field[f] = body_velocity(xc[0], xc[1], xc[2], c, s, poly, vel);
}

extern "C" __global__ void e3_geom_seam(GeomGrid g, const int* __restrict__ open_w, double* __restrict__ w)
{
    long long n = (long long) blockIdx.x * blockDim.x + threadIdx.x;
    long long nxy = (long long) g.nx * g.ny;
    if (n >= nxy) return;
    long long fn = (long long) g.nz * nxy + n;
    if (open_w[n] == 0 && open_w[fn] == 0) w[fn] = w[n];
}