Metadata-Version: 2.4
Name: gflex
Version: 2.0.0b1
Summary: One- and two-dimensional plate bending, designed for Earth's lithosphere
Author-email: "Andrew D. Wickert" <awickert@umn.edu>, "Eric W. H. Hutton" <eric.hutton@colorado.edu>
Maintainer-email: "Andrew D. Wickert" <awickert@umn.edu>
License: GPL-3.0-only
Project-URL: homepage, https://github.com/awickert/gflex
Project-URL: documentation, https://gflex.readthedocs.io
Project-URL: repository, https://github.com/awickert/gflex
Project-URL: changelog, https://github.com/awickert/gflex/blob/master/CHANGES.md
Keywords: geophysics,geology,geodynamics,lithosphere,isostasy,GRASS GIS
Classifier: Development Status :: 5 - Production/Stable
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: GNU General Public License v3 (GPLv3)
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: Implementation :: CPython
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.11
Description-Content-Type: text/markdown
License-File: LICENSE.md
License-File: AUTHORS.md
Requires-Dist: matplotlib
Requires-Dist: numpy
Requires-Dist: pyyaml
Requires-Dist: scipy
Provides-Extra: bmi
Requires-Dist: bmipy; extra == "bmi"
Provides-Extra: plot
Requires-Dist: cmcrameri; extra == "plot"
Provides-Extra: dev
Requires-Dist: nox; extra == "dev"
Provides-Extra: testing
Requires-Dist: pytest; extra == "testing"
Provides-Extra: docs
Requires-Dist: sphinx>=9.1; extra == "docs"
Requires-Dist: furo>=2025.12.19; extra == "docs"
Requires-Dist: myst-parser>=5.0; extra == "docs"
Dynamic: license-file

[![CSDMS Component][csdms_badge]][csdms_gflex] [![DOI][doi_badge]][doi_link] [![Test][test_badge]][test_workflow]

# gFlex

***Multiple methods to solve elastic plate flexure, designed for applications to Earth's lithosphere.***

These instructions are meant to take a user familiar with computers but new to (or a beginner with) Python through the basics of how to get gFlex to work. The Python scripting part towards the end should be pretty straightforward as well, insofar as information is provided on how to get and set the chosen values inside gFlex. *Please leave a message if you have trouble working with gFlex; your comments could assist both you and the more general improvement of this documentation.*

When you use gFlex, please cite:

**Wickert, A. D. (2016), [Open-source modular solutions for flexural isostasy: gFlex v1.0][paper_doi], *Geosci. Model Dev.*, *9*(3), 997–1017, doi:10.5194/gmd-9-997-2016.**

If you additionally want an up-to-date citation for the latest source-code release, please see that given in [CITATION.cff](CITATION.cff).

## Documentation

Full documentation, including a configuration file parameter reference, accuracy benchmarks, and the full API, is available at **[gflex.readthedocs.io](https://gflex.readthedocs.io)**.

## Installation

```bash
pip install gflex
```

gFlex requires **Python ≥ 3.11**. Dependencies (numpy, scipy, matplotlib, pyyaml) are installed automatically.

For a development install from source:

```bash
git clone https://github.com/awickert/gFlex.git
cd gFlex
pip install -e .
```

## Running

Once gFlex is installed, it is possible to run it in four ways:
 1. With a configuration file
 2. Within a Python script
 3. Within GRASS GIS
 4. As part of the Landlab Earth-surface modeling framework, with a CSDMS Basic Model Interface (BMI)

For options 1 and 2, there are pre-built methods that can be selected along the way to visualize results. These use Python's Matplotlib plotting library. For option 3, GRASS GIS is used for visualization. In option 4, output from Landlab can be visualized with Matplotlib.

#### With configuration file

A configuration file can be generated to run gFlex; see examples in the **input/** directory. To run gFlex using this file, one simply opens a terminal window and types:

```bash
gflex <path-to-configuration-file>
```

This can be run from any directory, as the installation of gFlex adds the program "gflex" to the system path.

Configuration files must be YAML (`.yaml` / `.yml` extension). See **input/input_f1d.yaml** and **input/input_f2d.yaml** for complete examples.

A minimal 2-D YAML configuration file looks like:

```yaml
mode:
  dimension: 2
  method: fd
parameter:
  youngs_modulus: 6.5e10
  poissons_ratio: 0.25
  gravitational_acceleration: 9.8
  mantle_density: 3300
  infill_material_density: 0
input:
  loads: path/to/loads.txt
  elastic_thickness: path/to/Te.txt
output:
  plot: both
numerical:
  grid_spacing_x: 4000
  boundary_condition_west: zero_moment_zero_shear
  boundary_condition_east: zero_displacement_zero_slope
numerical2D:
  grid_spacing_y: 4000
  boundary_condition_north: mirror
  boundary_condition_south: zero_slope_zero_shear
```

For a full parameter reference, see the [Configuration Files](https://gflex.readthedocs.io/en/latest/configuration.html) page on ReadTheDocs.

#### Finite-difference boundary conditions

Six boundary conditions are available for FD solutions (see also Table 1 in Wickert, 2016):

| Name | Condition | Physical interpretation |
|------|-----------|------------------------|
| `zero_displacement_zero_slope` | w = 0, dw/dx = 0 | Clamped end: zero deflection and zero slope (no rotation) |
| `zero_displacement_zero_moment` | w = 0, d²w/dx² = 0 | Simply supported (pinned) end: zero deflection, free to rotate |
| `zero_moment_zero_shear` | d²w/dx² = d³w/dx³ = 0 | Broken plate: free end with no moment or shear |
| `zero_slope_zero_shear` | dw/dx = d³w/dx³ = 0 | Plate is level at the boundary but free to deflect there; no shear transmitted |
| `mirror` | w(b − x) = w(b + x) | Mirror-symmetry plane — model only half of a symmetric system |
| `periodic` | w(0) = w(L) | Wrap-around: the domain tiles infinitely in both directions |

For SAS and SAS_NG, `no_outside_loads` (or a blank entry) is used instead; the plate is assumed undeflected at infinity.

**FD boundary-condition warnings:** when running F1D or F2D with the finite-difference solver, gFlex issues `UserWarning` messages for `'zero_moment_zero_shear'` (free broken plate end — verify a rifted margin is intended), `'zero_slope_zero_shear'` (no clear geological analog), and when the nearest loaded cell is within one flexural wavelength of a `'zero_displacement_zero_slope'` boundary (the forebulge would be suppressed). See the [API reference](https://gflex.readthedocs.io/en/latest/api.html#fd-boundary-condition-warnings) for how to suppress or re-enable these warnings.

**A note on `zero_slope_zero_shear`:** the label in Wickert (2016) is "free displacement of a horizontally clamped boundary." The plate is forced to be exactly level at the boundary (dw/dx = 0) — as if it were clamped against rotation — while its vertical position is unconstrained and no shear force is transmitted (d³w/dx³ = 0). This is superficially similar to `mirror` (both enforce zero slope at the boundary), but `zero_slope_zero_shear` uses a different finite-difference stencil and the two produce noticeably different solutions even far from the boundary. For symmetry problems — e.g., modelling half of a symmetric mountain range or ice sheet — `mirror` is the more accurate choice.

#### Within a Python script (with or without a configuration file)

You may run gFlex from other Python programs. When you install it (above), this also produces a Python module that you may import to access it while scripting.

##### With no configuration file (recommended)
**input/run_in_script_2D.py**, reproduced below, is a good example of how to set the variables and run the model. This method requires no input file, as all of the values are set inside the Python script that imports gflex. This is essentially how the GRASS GIS interface was written, and is a way to embed the abilities of gFlex into another model. A one-dimensional example, **input/run_in_script_1D.py**, is also available.

```python
#! /usr/bin/env python

import gflex
import numpy as np

flex = gflex.F2D()

flex.quiet = False

flex.method = 'fd' # Solution method: * fd (finite difference)
                   #                  * fft (spectral)
                   #                  * sas (superposition of analytical solutions)
                   #                  * sas_ng (ungridded SAS)

flex.g = 9.8 # acceleration due to gravity
flex.E = 65E9 # Young's Modulus
flex.nu = 0.25 # Poisson's Ratio
flex.rho_m = 3300. # mantle_density
flex.rho_fill = 0. # infill_material_density

flex.T_e = 35000.*np.ones((50, 50)) # Elastic thickness [m] -- scalar but may be an array
flex.T_e[:,-3:] = 0.
flex.qs = np.zeros((50, 50)) # Template array for surface load stresses
flex.qs[10:40, 10:40] += 1E6 # Populating this template
flex.dx = 5000. # grid cell size, x-oriented [m]
flex.dy = 5000. # grid cell size, y-oriented [m]
# Boundary conditions can be:
# (FD): zero_displacement_zero_slope, zero_displacement_zero_moment, zero_moment_zero_shear, zero_slope_zero_shear, mirror, or periodic
# For SAS or SAS_NG, no_outside_loads is valid, and no entry defaults to this
flex.bc_west = 'zero_displacement_zero_slope' # west boundary condition
flex.bc_east = 'zero_moment_zero_shear' # east boundary condition
flex.bc_south = 'zero_displacement_zero_slope' # south boundary condition
flex.bc_north = 'zero_displacement_zero_slope' # north boundary condition

# latitude/longitude solutions are exact for SAS, approximate otherwise
#latlon = # true/false: flag to enable lat/lon input. Defaults False.
#planetary_radius = # radius of planet [m], for lat/lon solutions

# Optional: in-plane stresses [Pa] (supported by FD and FFT solvers)
#flex.sigma_xx = 0.  # east–west compression/tension
#flex.sigma_yy = 0.  # north–south compression/tension
#flex.sigma_xy = 0.  # shear — couples x and y deflection

flex.initialize()
flex.run()

# Read deflection BEFORE calling finalize (finalize clears w and qs)
deflection = flex.w  # assign to a variable for later use

# If you want to plot or save output, do it before finalize:
flex.plot_choice='both'
# An output file for deflections could also be defined here
# flex.w_out_file =
flex.output() # Plots and/or saves output, or does nothing, depending on
              # whether flex.plot_choice and/or flex.w_out_file have been set

flex.finalize()
```

##### With a configuration file

If you would like to use a Python script with a configuration file, this is also possible.

```python
import gflex

# Pass the YAML file to the constructor; F1D for 1-D, F2D for 2-D problems.
filename = 'input/input_f1d.yaml'
flex = gflex.F1D(filename)

flex.initialize()
flex.run()

# Read deflection before finalize (finalize clears w and qs)
deflection = flex.w

# Standalone plotting output if you so desire
flex.plot_choice = 'w'
flex.output()

flex.finalize()
```


#### Within GRASS GIS

To run gFlex inside of GRASS GIS 8, install the addons from within a GRASS GIS session:

```bash
g.extension r.flexure
g.extension v.flexure
```

**r.flexure** is used for raster grids by either finite difference or analytical methods. **v.flexure** takes advantage of the ungridded analytical method to solve for flexure at an arbitrary set of load points, albeit limited to cases with constant elastic thickness. The source code and manual pages are available in the [GRASS GIS addons repository](https://github.com/OSGeo/grass-addons).

When running **r.flexure**, it is important to ensure that the elastic thickness map is at or properly interpolated to the computational region (**g.region**) resolution before solving. A nearest-neighbor interpolated Te map will cause perceived gradients in elastic thickness to be very sharp, and this will strongly affect (and misdirect) the flexural solutions.

#### As part of Landlab and CSDMS

[Landlab](https://landlab.github.io) is an Earth-surface modeling framework built to facilitate easy integration of geomorphic, ecological, hydrological, geological, and other Earth-surface models. gFlex can be used as a Landlab component; see the [Landlab repository](https://github.com/landlab/landlab) for details.

gFlex also implements the [CSDMS Basic Model Interface (BMI)](https://bmi.readthedocs.io), enabling it to be coupled with other models in the CSDMS framework. The BMI wrapper is available as `gflex.BmiGflex` and requires the optional `bmipy` dependency:

```bash
pip install gflex[bmi]
```

### Plotting

There are four plot choices, defined via `self.plot_choice`:
* `'q'`: plots the load in mantle-density-equivalent units of length
* `'w'`: plots the deflection in units of length
* `'both'`: plots both deflection and loads in separate panels of a 2-subplot figure
* `'combo'`: (1D only): plots lithospheric deflections and the deflected mantle-density-equivalent load atop it.
  * Note that the load does not affect the area above/below the datum filled when `rho_fill != 0`. This affects the buoyant balance associated with the motion of the plate, with no additional considerations for topography. If you would like to include topography, an iterative approach (e.g., finding areas below sea level, filling them, flexing, finding new areas below sea level, and so on) is recommended.

## Utilities

The **input/** directory contains several example scripts. The public API also provides standalone utility functions importable from `gflex`:

* `flexural_wavelengths(Te, ...)` — computes the flexural parameter α, first zero-crossing, and flexural wavelength for a given elastic thickness; useful for choosing grid spacing and domain size.

Domain-padding utilities (reduce spurious boundary effects when using variable *Te*):

*2-D (F2D):*
* `recommended_pad_width(Te, dx, ...)` — returns the recommended padding width (in cells).
* `smooth_pad_Te(Te, pad_width, ...)` — extends a 2-D variable-*Te* array with a smooth linear taper.
* `pad_domain(Te, qs, dx, ...)` — pads both the 2-D elastic thickness and load arrays and returns the padding width.

*1-D (F1D):*
* `recommended_pad_width_1d(Te, dx, ...)` — returns the recommended 1-D padding width (in cells).
* `smooth_pad_Te_1d(Te, pad_width, ...)` — extends a 1-D variable-*Te* array with a smooth linear taper.
* `pad_domain_1d(Te, qs, dx, ...)` — pads both the 1-D elastic thickness and load arrays and returns the padding width.

See the [API reference](https://gflex.readthedocs.io/en/latest/api.html) for full documentation of these functions.


[csdms_badge]: https://custom-icon-badges.demolab.com/badge/CSDMS-Component-2473c2?logo=csdms&style=for-the-badge
[csdms_gflex]: https://csdms.colorado.edu/wiki/Model:GFlex
[doi_badge]: https://zenodo.org/badge/DOI/10.5281/zenodo.10471939.svg
[doi_link]: https://doi.org/10.5281/zenodo.10471939
[test_badge]: https://github.com/awickert/gflex/actions/workflows/test.yml/badge.svg
[test_workflow]: https://github.com/awickert/gflex/actions/workflows/test.yml
[paper_doi]: https://www.geosci-model-dev.net/9/997/2016/gmd-9-997-2016.html

# Credits

## Development Lead

* [Andrew D. Wickert](https://github.com/awickert)

## Contributors

* [Eric W. H. Hutton](https://github.com/mcflugen)

# Release Notes

## [2.0.0b1] - 2026-06-04

### Breaking changes

- **`finalize()` now clears all model state** — ``w``, ``qs``, and the
  coefficient matrix are deleted.  Read ``w`` and call ``output()``
  **before** ``finalize()``.  The CLI (``gflex.py``) has been fixed to
  follow this order.

- **INI configuration file format removed** — only YAML (`.yaml` / `.yml`)
  configuration files are supported.  Passing a `.ini` file now raises
  `ValueError`.  Convert any remaining INI configs to YAML format.

- **Iterative FD solver removed** — `F1D` now always uses a direct sparse LU
  factorization.  The `ConvergenceTolerance` configuration key is ignored.
- **G2009 plate solution removed** — the finite-difference solver now always
  uses the vWC1994 stencil.  Setting `PlateSolutionType = 'G2009'` previously
  selected a different stencil; it is now a no-op attribute (ignored silently).
- **`PlateSolutionType` removed from the public interface** — the attribute is
  no longer documented or read by any solver.  Existing scripts that set it
  will not raise an error, but the value has no effect.
- **`te` attribute renamed to `T_e`** — the elastic-thickness attribute on
  `F1D` and `F2D` is now `T_e` (PEP 8 subscript notation; avoids ambiguity
  with a time variable `t` and a generic elastic parameter `e`).  `flex.te`
  now raises `AttributeError`.  Update all occurrences: `flex.te = …` → `flex.T_e = …`.

- **`"mirror"` is now an alias for `"zero_slope_zero_shear"`** — the canonical
  BC name is `"zero_slope_zero_shear"` (consistent with all other long-form names);
  `"mirror"` remains a permanent short alias (like `"clamped"` and `"free"`) and
  does not trigger any warning.

- **`solver` attribute now raises `ValueError` for unsupported values** —
  only `"direct"` is supported in this release.  Passing any other string
  (e.g. `"iterative"`) previously silently fell back to the direct solver;
  it now raises `ValueError` immediately.  An iterative solver may be added
  in a future version, at which point additional values will become valid.

- **BC string case normalised; v1.x PascalCase strings removed** — all
  boundary-condition strings are now lowercase (`"zero_displacement_zero_slope"`,
  `"zero_moment_zero_shear"`, `"zero_slope_zero_shear"`, `"periodic"`).  The old v1.x
  PascalCase names (`"0Displacement0Slope"`, `"0Moment0Shear"`, `"Mirror"`,
  etc.) now raise `ValueError`.  See the `boundary_conditions` page for the
  full mapping.

### New features

- **Inhomogeneous (prescribed-value) boundary conditions** in 1-D and 2-D FD
  solver.  Any edge can carry a combination of prescribed displacement, slope,
  moment, and/or shear by passing a dict instead of a string:

  ```python
  flex.bc_west = {"displacement": w_arr, "slope": dwdx_arr}
  flex.bc_east = {"moment": M0, "shear": V0}
  ```

  The primary use case is **nested-domain (coarse-to-fine) modelling**: extract
  `w` and `np.gradient`-based slopes at the boundary of a fine sub-domain from
  a coarse regional solve, impose them as BCs, and the fine solver inherits the
  full far-field loading through those four boundary arrays.  Also supports
  broken-plate problems with applied edge loads.  Both 1-D and 2-D; verified
  by a round-trip unit test (`TestNestedModelGradientRoundTrip`).

- **`gflex.VALID_BC_STRINGS_1D` and `gflex.VALID_BC_STRINGS_2D`** — public
  `frozenset` constants listing every accepted BC string (canonical names and
  aliases).  Wrappers can validate user input against these sets instead of
  maintaining a parallel copy that drifts with new releases.

- **`"clamped"`, `"free"`, and `"mirror"` BC aliases** — concise alternatives to the full
  canonical names: `"clamped"` normalises to `"zero_displacement_zero_slope"`;
  `"free"` normalises to `"zero_moment_zero_shear"`;
  `"mirror"` normalises to `"zero_slope_zero_shear"`.  All three produce
  bit-identical results to their canonical names.

### Performance

- **LU factorization cache** — set ``flex.cache_factorization = True`` to
  cache the sparse-LU factorization of the FD coefficient matrix across
  ``run()`` calls.  The cache is reused automatically when the matrix is
  unchanged (Te, grid, BCs, and physical parameters all fixed) and
  invalidated by a hash comparison when any of those inputs change.  Set
  ``flex.cache_factorization = "no_check"`` to skip the hash and reuse the
  factorization unconditionally — maximum performance for coupling workflows
  where the user guarantees matrix stability.  Default ``False`` preserves
  existing behaviour.
- **Multithreaded FFT** — ``scipy.fft`` transforms now use all available CPU
  threads (``workers=-1``).  Benefit is marginal at typical grid sizes;
  meaningful at very large grids (≳ 2000 × 2000 cells).

### Improvements

- `twoSurfplots()` now shows three panels (load, elastic thickness, deflection)
  when `Te` is a 2-D array, using an asymmetric diverging colormap for the
  deflection.
- `export_for_blender()` writes load-cylinder geometry (radius from loaded
  area, height proportional to load magnitude) and per-vertex colour weights
  to the mesh file, enabling the companion Blender script to render a
  physically scaled load cylinder.
- Blender scene script: cylinder base deforms with the plate surface; camera
  frames the full scene vertically.

### Bug fixes

- Fixed `SyntaxWarning` for invalid escape sequences (`\sigma`) in the F2D
  class docstring (Python 3.12+).
- Fixed a longstanding ghost-node error in the 1-D
  `zero_displacement_zero_slope` FD boundary condition.  The original
  implementation (present since the first version) dropped ghost-node stencil
  terms rather than eliminating them via even reflection, and left the
  boundary row coupled to interior nodes rather than decoupling it as a
  Dirichlet constraint.  The practical effect was that the boundary deflection
  was not strictly zero and the zero-slope condition was not enforced; results
  were accurate only when the boundary was far enough from any load that
  deflection was naturally negligible there (the recommended use case, and the
  intent of the existing proximity warning).  The corrected implementation
  decouples the boundary row (enforcing w = 0 exactly) and folds the
  even-reflected ghost into the adjacent interior node (enforcing dw/dx = 0).
- Fixed the same two ghost-node bugs in the **2-D `zero_displacement_zero_slope`**
  boundary condition (boundary rows/columns not decoupled as Dirichlet
  constraints; even-reflection ghost absent at the first interior row/column on
  all four edges).  Convergence recovers from first order (O(Δx^0.92)) to
  second order (O(Δx^1.99)).
- Fixed a ghost-node inconsistency in the **`zero_moment_zero_shear`** FD
  boundary condition (1-D and 2-D).  The first interior node's stencil used a
  shear ghost evaluated at a staggered location one cell inward, inconsistent
  with the moment ghost used at the boundary node itself.  The corrected
  implementation uses the moment condition consistently at both rows, recovering
  second-order convergence (O(Δx^2.00) in 1-D, O(Δx^2.01) in 2-D) from first
  order.

### Documentation

- New `boundary_conditions` page: comprehensive reference covering all BC
  types with a summary table (structural-mechanics and geophysical names), SVG
  diagrams for each condition, physical guidance on when each is appropriate,
  and a deprecation table mapping v1.x PascalCase strings to their v2.0
  lowercase equivalents.
- New `greenland_example` page: a realistic Greenland ice-sheet flexure model
  (BedMachine v6 ice thickness; Steffen et al. spatially variable elastic
  thickness) and a hypothetical nested-domain seamount scenario that illustrates
  how inhomogeneous BCs couple a coarse regional model to a fine local one.
- Accuracy page extended with MMS convergence tables and figures for the 2-D
  `zero_displacement_zero_slope` ghost-node fix and the 1-D/2-D
  `zero_moment_zero_shear` ghost-node fix.

### Tests

- 335 tests passing across 1-D and 2-D FD, FFT, SAS/SAS_NG solvers, all BC
  types, inhomogeneous BCs, domain padding, warnings, and BC aliases.
- `TestNestedModelGradientRoundTrip`: verifies that `np.gradient`-extracted
  slopes used as inhomogeneous Dirichlet BCs reproduce the full-domain interior
  to within 2 % of peak deflection, confirming the slope sign convention for
  all four edges in the nested-domain workflow.
- Parametrised alias tests confirm `"clamped"` and `"free"` produce
  bit-identical results to their canonical names in both 1-D and 2-D.

## [1.4.0](https://github.com/awickert/gFlex/releases/tag/v1.4.0) - 2026-05-29

### New features

- **FFT spectral solver** for 1-D and 2-D problems.  Exact for periodic
  boundaries; other boundary conditions are handled with 4α zero-padding
  (NoOutsideLoads approximation).  Requires scalar (uniform) elastic
  thickness.  In-plane stresses (`sigma_xx`, `sigma_yy`, `sigma_xy`) are
  fully supported.
- **YAML configuration file format** alongside the legacy INI format.  Pass
  a `.yaml` / `.yml` file to the CLI or to the `F1D` / `F2D` constructor.
- **Domain-padding utilities** to reduce spurious boundary effects when
  using spatially variable elastic thickness:
  - 2-D: `pad_domain`, `smooth_pad_Te`, `recommended_pad_width`
  - 1-D: `pad_domain_1d`, `smooth_pad_Te_1d`, `recommended_pad_width_1d`
- **FD boundary-condition warnings** — `F1D` and `F2D` now issue
  `UserWarning` messages for `'0Moment0Shear'` (free broken end — verify
  a rifted margin is intended), `'0Slope0Shear'` (no clear geological
  analog), and when the nearest loaded cell is within one flexural
  wavelength of a `'0Displacement0Slope'` boundary (forebulge suppression).

### Improvements

- Iterative FD solver (`Solver = iterative`) upgraded from Jacobi to ILU
  preconditioning; falls back to a direct solve if LGMRES does not converge.
  The `ConvergenceTolerance` parameter is now passed as `rtol` (relative
  residual) to SciPy's LGMRES; default remains `1e-3`.
- In-plane membrane stresses (`sigma_xx`, `sigma_yy`, `sigma_xy`) now
  supported by all FD and FFT solvers in both 1-D and 2-D.

### Bug fixes

- Fixed 2-D FD constant-Te stencil: corrected dx/dy swap, σ swap, and
  missing `/dx²` factors that produced incorrect results for asymmetric grids.
- Fixed 2-D FFT `sigma_xy` assembly: double-wrapped corner entries and wrong
  coefficient in the Periodic right-roll buffer.
- Fixed `0Slope0Shear` description in configuration reference and README.

### Documentation

- New Sphinx / Read the Docs site with a theory page (governing equations,
  solution methods, in-plane stresses), an accuracy page, a configuration
  reference, API reference, and changelog.
- Theory page covers the full governing PDE with in-plane stress terms for
  both 1-D and 2-D, the variable-D FD expansion, and the FFT transfer
  function.

### Tests

- Comprehensive new test suites for 1-D SAS/SAS_NG, 2-D SAS/SAS_NG,
  1-D FFT, 2-D FFT (including sigma_xy), 1-D FD boundary conditions with
  analytical cross-validation, 2-D FD, FD boundary-condition warnings, and
  all domain-padding utilities.

## [1.3.0](https://github.com/awickert/gFlex/releases/tag/v1.3.0) - 2026-05-27

- New: `BmiGflex` — CSDMS Basic Model Interface (BMI) v2 implementation;
  `bmipy` is an optional dependency and gFlex works without it.
- New: `F1D`, `F2D`, and `BmiGflex` are now accessible directly from the
  `gflex` package namespace (e.g. `import gflex; gflex.F2D()`); previously
  this raised `AttributeError`.
- Fix: updated `scipy.sparse.linalg` import for compatibility with modern
  scipy.
- Fix: latent typo `self.T_e` → `self.Te` in the scalar-Te branch of
  `get_coeff_values` (introduced 2021-06-26); would produce incorrect
  results when tectonic stress terms (sigma_xx, sigma_yy, sigma_xy) are
  non-zero.
- Removed outdated root-level `gflex_bmi.py`.

## [1.2.0](https://github.com/awickert/gFlex/releases/tag/v1.2.0) - 2024-01-08

- Python 3 modernisation: black formatting, isort, updated imports,
  argparse CLI, `pyproject.toml` replacing `setup.py`.
- GitHub Actions CI replacing Travis CI.
- Eric W. H. Hutton added as author.

## [1.1.1](https://github.com/awickert/gFlex/releases/tag/v1.1.1) - 2021-06-26

- Updated PyPI support: twine upload, README.md
- Throw meaningful error if a nonuniform self.Te grid is used with the
  analytical solution

## [1.1.0](https://github.com/awickert/gFlex/releases/tag/v1.1.0) - 2018-05-28

- Support for both Python 2 and Python 3
- Main code
- Examples
- README.md updates
- PATH updates in
- Code tests included
- Code testing on commit by Travis
- Updated on PyPI

## [1.0.1](https://github.com/awickert/gFlex/releases/tag/v1.0.1) - 2018-05-25

- Final Python 2 (only) release
- Additional documentation added

## [1.0.0](https://github.com/awickert/gFlex/releases/tag/v1.0.0) - 2017-05-10

This is the update is what is available now from pypi.

- Minor updates to the main gFlex codes
- Addition of a missing "12" in the flexural wavelength calculator in the utilities.


## [1.0](https://github.com/awickert/gFlex/releases/tag/v1.0) - 2016-01-28

- First full release of gFlex in association with the now-accepted GMD paper,
"Open-source modular solutions for flexural isostasy: gFlex v1.0", by A. D. Wickert.


## [0.9](https://github.com/awickert/gFlex/releases/tag/v0.9) - 2015-03-05

- Release submitted to GMD and that appears on PyPI; this is the same as v1.0a
  on umn-earth-surface.

## [0.8.1](https://github.com/awickert/gFlex/releases/tag/v0.8.1) - 2015-03-04

- Fixed error in PyPI integration from v0.8 and updated README.md to include
  PyPI integration.


## [0.8](https://github.com/awickert/gFlex/releases/tag/v0.8) - 2015-03-04

- First fully-functional and error-checked release, and therefore the first
  true non-"pre-release".


## [0.7](https://github.com/awickert/gFlex/releases/tag/0.7) - 2015-01-21

This release may be short-lived, but is the turning point at which work
on gFlex will go from fixing deficiencies to adding new capabilities.
As such, it is receiving version 0.7 (the last version of its parent project,
"Flexure", was 0.6 back in 2012). It now seems to be functional and stable,
and thus this tagged release will be the basis for this stage of the model
to be pulled to @csdms-contrib and @umn-earth-surface.
