================================================================================
 final_TPMS - Triply Periodic Minimal Surface (TPMS) Porous Geometry Generator
================================================================================

OVERVIEW
--------
This MATLAB code builds triangulated surface meshes of TPMS-based porous
media (Schwarz P, Gyroid, Schwarz D and variants) at a chosen porosity,
tiles the unit cell into a block, and can export the result as an STL file
for meshing, CFD or 3D printing. It can also compute surface functionals of
the unit cell (surface area, area-averaged mean and Gaussian curvature, and
cap areas on the cell faces), and plot the curvature on the surface.

A separate script, pressure_raw_data.m, stores measured or digitised
pressure-gradient vs velocity data for several of these geometries.


HOW IT WORKS
------------
  1. The TPMS is defined as an implicit level set  u(x,y,z) = rhs  on a
     unit cube sampled at geometry.resolution points per edge.
  2. MATLAB's isosurface extracts the surface (the "iso" part) and isocaps
     closes it on the cube faces (the "caps"), giving a closed solid.
  3. The unit cell is scaled by Lref / anisotropy in each direction.
  4. getCells repeats the cell Lxyz/Lref times in x, y and z. The caps are
     only kept on the outer faces of the block (copied from the min-x,
     min-y and min-z caps, which is valid because the surface is periodic).
  5. The surface and caps are merged into one face/vertex set (fTot, vTot)
     and optionally written to STL.


SUPPORTED GEOMETRIES
--------------------
  geometry.name   Level-set function u                     Porosities
  -------------   ---------------------------------------  ------------------
  "P"             Schwarz Primitive:                       0.25, 0.5, 0.65,
                  cos X + cos Y + cos Z                    0.75
  "Pi"            Inverted P (same u, opposite level and   0.5, 0.65, 0.75
                  caps on the other side, so solid and
                  void are swapped)
  "G"             Gyroid:                                  0.5, 0.65, 0.75
                  sinX cosY + sinY cosZ + sinZ cosX
  "D"             Schwarz Diamond (4-term form)            0.5, 0.65, 0.75
  "D_alt"         Schwarz Diamond (2-term form):           0.5, 0.65, 0.75
                  cosX cosY cosZ - sinX sinY sinZ

  (X = 2*pi*x, etc.) Each porosity maps to a hard-coded level-set constant
  rhs in getSurfacesTPMS.m. Any other name/porosity combination is not
  supported (rhs is left undefined and MATLAB will error).


FILES
-----
TPMS_generator.m
    Driver script. Set the geometry and the operation flags at the top,
    then run it. This is the file to run.

getSurfacesTPMS.m
    [faces, vertices] = getSurfacesTPMS(geometry)
    Builds one unit cell. Returns structs with:
      .iso   - faces/vertices of the TPMS surface
      .cap   - struct with .xmin, .ymin, .zmin caps (used for tiling)
      .caps  - caps on all six cube faces (used for the functionals and
               the single-cell plot)

getCells.m
    [fTot, vTot] = getCells(faces, vertices, nbr, Lref, cap)
    Tiles a face/vertex set into an nbr(1) x nbr(2) x nbr(3) block.
    With cap = false it copies the surface; with cap = true it builds the
    caps on the six outer faces of the block from the .xmin/.ymin/.zmin
    caps.

initPlots.m
    Sets MATLAB graphics defaults (LaTeX interpreter, Times New Roman,
    font size 12, line width 1).

patchcurvature.m
    Third-party function by D. Kroon (University of Twente, 2011; last
    updated 2014), originally from the MATLAB File Exchange. Computes
    per-vertex mean and Gaussian curvature and principal directions of a
    triangulated mesh by local quadratic fitting. Self-contained.

modstlwrite.m
    Modified version of stlwrite by Sven Holcombe (2011, MATLAB File
    Exchange), last modified by Thomas Hunter. Writes binary (default) or
    ASCII STL from faces/vertices, an FV struct, or gridded X/Y/Z data.
    The target folder must already exist.

pressure_raw_data.m
    Standalone data script (not called by the others). Defines one struct
    per geometry - P (P75), P65, P50, D (D75), D65, D50, G (G75) and
    Pi (Pi75) - each with:
      .L              a length [m] (2 mm or 6.35 mm)
      .vt,  .dpdt     one velocity / pressure-gradient series
      .vT,  .dpdT     a second velocity / pressure-gradient series
    t and T refer to the medium thickness (respectively, t=2L and T=10L).
    The data is in SI units.

REQUIREMENTS
------------
  - MATLAB R2022a or newer (clim is used in the curvature plots; on older
    versions replace clim with caxis).
  - No additional toolboxes.


HOW TO RUN
----------
  1. Put all .m files in the same folder (or on the MATLAB path).
  2. Edit the settings at the top of TPMS_generator.m.
  3. Run:
         >> TPMS_generator
  4. The final geometry is left in the workspace as fTot (faces) and
     vTot (vertices).

  To export an STL, create a folder called stl_files next to the scripts,
  set save_surface = true, and fix the file-name bug described under
  Known issues (item 1).


SETTINGS IN TPMS_generator.m
----------------------------
  geometry.name         "Pi"          TPMS type (see Supported geometries)
  geometry.Lref         1 m           Unit-cell edge length
  geometry.Lxyz         [1 1 1]*Lref  Size of the block; Lxyz/Lref must be
                                      whole numbers (cells per direction)
  geometry.resolution   64            Grid points per unit-cell edge used
                                      to extract the isosurface
  geometry.anisotropy   [1 1 1]       Divides the cell size in each
                                      direction (see Known issues, item 3)
  geometry.porosity     0.75          Target porosity (see table above)

  plot_surface          false         Plot the single cell (surface + caps)
  save_surface          false         Write the tiled geometry to STL
  save_name             "temp"        Intended STL file name
  compute_functionals   false         Compute the surface functionals below
  plot_functionals      false         Plot mean and Gaussian curvature on
                                      the surface (needs
                                      compute_functionals = true)


SURFACE FUNCTIONALS (compute_functionals = true)
------------------------------------------------
Computed for a single unit cell:
  M_0      porosity (taken directly from the input)
  M_1      total area of the TPMS surface
  M_2      area-weighted mean of the mean curvature
  M_3      area-weighted mean of the Gaussian curvature
  barA_i   total cap area on the six cell faces
  A_i      total cell-face area minus the cap area

Per-face curvatures are the average of the three vertex values from
patchcurvature (with third-order neighbours), and only the real part is
used.


KNOWN ISSUES / THINGS TO CHECK
------------------------------
  1. Vertices are not merged where tiled cells meet, so the STL has
     duplicate vertices along the seams. Some meshing tools may need the
     mesh cleaned or merged first.