# MAC Grid Conventions

This document describes the coordinate and indexing conventions used throughout the PolyCFD codebase for the Marker-And-Cell (MAC) staggered grid.

## Overview

PolyCFD uses a standard MAC (Marker-And-Cell) staggered grid where:
- **Pressure** (and other scalar fields) are stored at **cell centers**
- **Velocity components** are stored at **face centers**, staggered in their respective directions

## Coordinate System

The domain spans from `(0, 0, 0)` to `(Lx, Ly, Lz)` where:
- `Lx = Nx × dx`
- `Ly = Ny × dy`  
- `Lz = Nz × dz`

Grid indices range:
- Cell indices: `i ∈ [0, Nx-1]`, `j ∈ [0, Ny-1]`, `k ∈ [0, Nz-1]`
- U-face indices: `i ∈ [0, Nx]`, `j ∈ [0, Ny-1]`, `k ∈ [0, Nz-1]`
- V-face indices: `i ∈ [0, Nx-1]`, `j ∈ [0, Ny]`, `k ∈ [0, Nz-1]`
- W-face indices: `i ∈ [0, Nx-1]`, `j ∈ [0, Ny-1]`, `k ∈ [0, Nz]`

## Position Conventions

### Cell Centers (Pressure, Scalars)

Cell `(i, j, k)` has its center at:

```
x = (i + 0.5) × dx
y = (j + 0.5) × dy
z = (k + 0.5) × dz
```

### Velocity Face Centers

#### U-Velocity (X-component)

U-face `(i, j, k)` is located at:

```
x = i × dx                    // Face normal to X, at integer i
y = (j + 0.5) × dy           // Cell-centered in Y
z = (k + 0.5) × dz           // Cell-centered in Z
```

The U-face `(i, j, k)` lies on the **left** face of cell `(i, j, k)`, between cells `(i-1, j, k)` and `(i, j, k)`.

#### V-Velocity (Y-component)

V-face `(i, j, k)` is located at:

```
x = (i + 0.5) × dx           // Cell-centered in X
y = j × dy                    // Face normal to Y, at integer j
z = (k + 0.5) × dz           // Cell-centered in Z
```

The V-face `(i, j, k)` lies on the **bottom** face of cell `(i, j, k)`, between cells `(i, j-1, k)` and `(i, j, k)`.

#### W-Velocity (Z-component)

W-face `(i, j, k)` is located at:

```
x = (i + 0.5) × dx           // Cell-centered in X
y = (j + 0.5) × dy           // Cell-centered in Y
z = k × dz                    // Face normal to Z, at integer k
```

The W-face `(i, j, k)` lies on the **back** face of cell `(i, j, k)`, between cells `(i, j, k-1)` and `(i, j, k)`.

## Visual Diagram (2D Slice in X-Y Plane)

```
        j+1  ──────V(i,j+1)──────
              │                 │
              │                 │
              │    P(i,j)       │
        j     U(i,j)   ×    U(i+1,j)
              │                 │
              │                 │
              │                 │
        j-1  ──────V(i,j)───────
              i               i+1

Legend:
  ×   = Cell center (pressure P)
  U() = U-velocity face (horizontal arrows →)
  V() = V-velocity face (vertical arrows ↑)
```

## Array Indexing

### Dense Interior Indexing (No Halos)

All interior arrays use dense 3D indexing with stride order `(i, j, k)`:

```csharp
// Cell-centered: Nx × Ny × Nz
int cellIdx = i + Nx * (j + Ny * k);

// U-faces: (Nx+1) × Ny × Nz
int uIdx = i + (Nx + 1) * (j + Ny * k);

// V-faces: Nx × (Ny+1) × Nz
int vIdx = i + Nx * (j + (Ny + 1) * k);

// W-faces: Nx × Ny × (Nz+1)
int wIdx = i + Nx * (j + Ny * k);
```

### Array Sizes

| Array Type | Dimensions | Total Size |
|------------|------------|------------|
| Cell (P, scalars) | Nx × Ny × Nz | Nx·Ny·Nz |
| U-faces | (Nx+1) × Ny × Nz | (Nx+1)·Ny·Nz |
| V-faces | Nx × (Ny+1) × Nz | Nx·(Ny+1)·Nz |
| W-faces | Nx × Ny × (Nz+1) | Nx·Ny·(Nz+1) |

## Divergence and Gradient Operators

### Divergence

For cell `(i, j, k)`, the discrete divergence is:

```
div(u) = (flux_x+ - flux_x- + flux_y+ - flux_y- + flux_z+ - flux_z-) / V

where:
  flux_x+ = αx[i+1,j,k] × Ax[i+1,j,k] × U[i+1,j,k]
  flux_x- = αx[i,j,k]   × Ax[i,j,k]   × U[i,j,k]
  (similar for Y and Z)
  
  V = cell volume (accounting for cut-cell geometry)
  αx = face aperture (0 = blocked, 1 = fully open)
  Ax = face area
```

### Pressure Gradient

The pressure gradient at faces uses central differences:

```
∂p/∂x at U-face (i,j,k) = (P[i,j,k] - P[i-1,j,k]) / dx
∂p/∂y at V-face (i,j,k) = (P[i,j,k] - P[i,j-1,k]) / dy
∂p/∂z at W-face (i,j,k) = (P[i,j,k] - P[i,j,k-1]) / dz
```

## Cut-Cell Geometry

For embedded boundaries, each cell and face has associated geometric quantities:

### Cell Quantities
- **CellClass**: `Fluid`, `Solid`, or `Cut`
- **Volume**: Fluid volume within the cell (0 for Solid, reduced for Cut)
- **VolumeFraction**: Volume / (dx × dy × dz)

### Face Quantities
- **FaceClass**: `Open`, `Closed`, or `Cut`
- **Alpha (α)**: Aperture fraction (0 = fully blocked, 1 = fully open)
- **Area**: Effective open area for flux calculation
- **Normal**: Unit normal vector at cut faces (for free-slip BC)

### Face Classification Rules

| FaceClass | Alpha | Condition |
|-----------|-------|-----------|
| Open | 1.0 | Both adjacent cells are fluid |
| Closed | 0.0 | Face straddles solid-fluid interface OR inside solid |
| Cut | 0 < α < 1 | Face partially intersects geometry |

**Important**: When a face straddles a solid-fluid interface (one adjacent cell is Solid, the other is Fluid/Cut), the face is forced to `Closed` with `α = 0` regardless of the SDF-based classification. This ensures no flux through solid boundaries.

## Moving Wall Velocity

For moving geometry (rotation, translation), wall velocity is applied at `Closed` faces only:

```
v_wall = ω × r   (for rotation)
v_wall = v_trans (for translation)

where:
  ω = angular velocity vector
  r = position - rotation_center
```

Only `Closed` faces receive wall velocity because:
- At `Closed` faces, α = 0, so `flux = α × A × v = 0` regardless of v
- Setting velocity at `Cut` faces would create artificial divergence

## Code Examples

### Computing Wall Velocity at a U-Face

```csharp
// U-face (i,j,k) position
float x = i * dx;
float y = (j + 0.5f) * dy;
float z = (k + 0.5f) * dz;

// For rotation about center with angular velocity omega
Float3 r = new Float3(x - center.X, y - center.Y, z - center.Z);
float wallVelX = omega.Y * r.Z - omega.Z * r.Y;  // (ω × r)_x
```

### Interpolating U-Velocity at Arbitrary Position

```csharp
// Convert world position to grid coordinates
// U-face at (i*dx, (j+0.5)*dy, (k+0.5)*dz)
float iFloat = x / dx;           // No offset - U is at integer i
float jFloat = y / dy - 0.5f;    // Offset by 0.5 - U is cell-centered in j
float kFloat = z / dz - 0.5f;    // Offset by 0.5 - U is cell-centered in k

// Then use trilinear interpolation...
```

## References

1. Harlow & Welch (1965) - Original MAC method
2. Chorin (1968) - Projection method for incompressible flow
3. Fedkiw et al. (2001) - Cut-cell methods for embedded boundaries