Metadata-Version: 2.4
Name: scikit-fem-native
Version: 0.3.5
Summary: Native finite-element assembly for scikit-fem-style Python workflows
Keywords: scikit-fem,finite-element-method,finite-elements,assembly,scientific-computing,computational-mechanics,sparse-matrix
License-Expression: Apache-2.0
License-File: LICENSE
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Operating System :: MacOS
Classifier: Operating System :: Microsoft :: Windows
Classifier: Operating System :: POSIX :: Linux
Classifier: Programming Language :: C++
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3 :: Only
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Programming Language :: Python :: 3.13
Classifier: Programming Language :: Python :: 3.14
Classifier: Topic :: Scientific/Engineering
Classifier: Topic :: Scientific/Engineering :: Mathematics
Project-URL: Homepage, https://github.com/kevin-tofu/scikit-fem-native
Project-URL: Documentation, https://scikit-fem-native.readthedocs.io/en/latest/
Project-URL: Repository, https://github.com/kevin-tofu/scikit-fem-native
Project-URL: Issues, https://github.com/kevin-tofu/scikit-fem-native/issues
Requires-Python: >=3.10
Requires-Dist: numpy>=1.24
Requires-Dist: scipy>=1.10
Requires-Dist: meshio<6.0.0,>=5.3.5
Provides-Extra: io
Requires-Dist: meshio>=5.3; extra == "io"
Provides-Extra: test
Requires-Dist: pytest>=8; extra == "test"
Requires-Dist: scikit-fem>=10; extra == "test"
Requires-Dist: shapely>=2; extra == "test"
Requires-Dist: matplotlib>=3.8; extra == "test"
Requires-Dist: meshio>=5.3; extra == "test"
Provides-Extra: benchmark
Requires-Dist: scikit-fem>=10; extra == "benchmark"
Requires-Dist: shapely>=2; extra == "benchmark"
Requires-Dist: matplotlib>=3.8; extra == "benchmark"
Provides-Extra: docs
Requires-Dist: sphinx>=8; extra == "docs"
Requires-Dist: furo>=2024.8.6; extra == "docs"
Requires-Dist: myst-parser>=4; extra == "docs"
Requires-Dist: sphinx-copybutton>=0.5.2; extra == "docs"
Requires-Dist: sphinx-design>=0.6.1; extra == "docs"
Description-Content-Type: text/markdown

# scikit-fem-native: native assembly for scikit-fem-style Python workflows

[![PyPI](https://img.shields.io/pypi/v/scikit-fem-native.svg)](https://pypi.org/project/scikit-fem-native/)
[![Python](https://img.shields.io/pypi/pyversions/scikit-fem-native.svg)](https://pypi.org/project/scikit-fem-native/)
[![CI](https://github.com/kevin-tofu/scikit-fem-native/actions/workflows/ci.yml/badge.svg)](https://github.com/kevin-tofu/scikit-fem-native/actions/workflows/ci.yml)
[![Documentation Status](https://readthedocs.org/projects/scikit-fem-native/badge/?version=latest)](https://scikit-fem-native.readthedocs.io/en/latest/)
[![License: LGPL-3.0](https://img.shields.io/badge/license-LGPL--3.0-blue.svg)](https://github.com/kevin-tofu/scikit-fem-native/blob/main/LICENSE)

**scikit-fem-native** is an independent native finite-element assembly engine for
[scikit-fem](https://github.com/kinnala/scikit-fem)-style Python workflows.  It
accelerates reusable assembly, geometry, quadrature, and sparse-scatter kernels
while preserving readable Python weak forms and, for its supported feature
subset, scikit-fem-compatible DOF ordering and numerical results.

> [!IMPORTANT]
> The distribution is named `scikit-fem-native` and imported as `skfemntv`.  It is
> an independent project, not an official scikit-fem distribution or backend.

## Motivation

The project exists to improve the performance of Python assembly workflows
represented by scikit-fem without replacing their Python programming model.
Weak forms, nonlinear constitutive updates, contact formulations, and solver
policy remain ordinary Python components; stable and reusable numerical kernels
are moved to native code where doing so provides a practical benefit.

This boundary is becoming more important in the era of AI-assisted coding.
Application-specific nonlinear formulations can now be developed, reviewed,
and revised rapidly when they remain small, explicit Python components.
Libraries such as scikit-fem therefore become more valuable, not less: they
provide a transparent mathematical vocabulary while avoiding the cost of
embedding every new formulation in a monolithic compiled solver.

`scikit-fem-native` aims to combine that flexibility with faster assembly.  It is a
selectable backend for Python finite-element applications, not a separate owner
of their material models or physical formulations.

## Built on the work of scikit-fem

This project exists because
[scikit-fem](https://github.com/kinnala/scikit-fem) demonstrated how effective
a finite-element library can be when mathematical forms remain concise,
readable Python.  scikit-fem is implemented in pure Python, which makes its
environment straightforward to install, inspect, and reproduce across
development machines, notebooks, and CI systems.  Its small and composable API
makes the connection between a weak formulation and the assembled system
unusually clear.  That transparency is valuable for teaching and research, but
also for production engineering: new elements, multiphysics terms, nonlinear
material models, and contact formulations can be inspected and changed without
hiding the governing equations behind a large generated or compiled framework.
See the [scikit-fem documentation](https://scikit-fem.readthedocs.io/) for its
installation guide, examples, and API reference.

scikit-fem also provides the practical foundation on which `scikit-fem-native` is
built: its public abstractions, element conventions, DOF ordering, examples,
and extensive body of implementation work define the compatibility target of
this project.  `scikit-fem-native` does not seek to replace that contribution.  It
explores how selected assembly kernels can be accelerated while preserving the
Python-first model that makes scikit-fem useful.

We are grateful to the scikit-fem authors and contributors.  The usefulness of
this backend is a direct consequence of the careful design and sustained work
of that community.

## Installation

```bash
python -m pip install scikit-fem-native
```

The distribution is named `scikit-fem-native`; import it as `skfemntv`.

## Quick start

```python
import skfemntv
from skfemntv.helpers import dot, grad

mesh = skfemntv.MeshTet()
basis = skfemntv.Basis(mesh, skfemntv.ElementTetP1())

@skfemntv.BilinearForm
def diffusion(u, v, w):
    return dot(grad(u), grad(v))

matrix = skfemntv.asm(diffusion, basis)
```

Applications can select `skfem.asm` or `skfemntv.asm` at their backend
boundary while preserving the weak-form call site.  `skfemntv` is not yet a
complete replacement for every scikit-fem feature, so backend selection should
remain explicit.

## Experimental H(div): affine RT0 simplices

Affine triangle and tetrahedron RT0 bases are available through the dedicated
`AffineTriRT0Basis` / `AffineTetRT0Basis` API. `SimplexRT0Assembler` assembles
vector mass, divergence-divergence, and vector-load terms. Cell quadrature is
computed by the C++ backend; sparse facet-DOF accumulation currently uses
Python/SciPy. This experimental path is separate from `Basis` and `asm`.

```python
import skfemntv

mesh = skfemntv.MeshTri.init_tensor([0., 0.5, 1.], [0., 0.5, 1.])
basis = skfemntv.AffineTriRT0Basis(mesh, intorder=4)
assembler = skfemntv.SimplexRT0Assembler(basis)
mass = assembler.assemble(value=1.0)
divergence = assembler.assemble(divergence=1.0)
```

## Example analysis: tension of a plate with a hole

This linear plane-stress example fixes horizontal motion at the left edge,
prescribes a small horizontal displacement at the right edge, and leaves
transverse motion free apart from one pinned DOF.  The image shows von Mises
stress normalized by the nominal stress.  With this mesh, the largest
cell-averaged stress is about 3.26 times nominal stress; the plotted value is
an element average, not a pointwise peak.  Refining the angular/radial mesh
from `180×20` to `360×40`, `720×80`, and `1080×120` gave peak ratios of
`3.231`, `3.263`, `3.275`, and `3.278`, respectively.  The reference value
`3` is for an infinite plate with a circular hole; this finite plate has a
hole diameter equal to 36% of its height, so its boundaries affect the peak.

![Von Mises stress around the hole in a tensioned plate](examples/plate_with_hole_tension.png)

## Example analysis: electrostatic potential

The same finite-element workflow also applies to a different physical model:
electrostatics.  Structural mechanics solves for a vector displacement field
`u` using an elasticity operator.  Electrostatics solves for a scalar potential
`phi` using the Poisson/Laplace operator
`-div(epsilon grad(phi)) = rho`, with electrical potential or charge-flux
boundary conditions.  The governing equation and unknown are different, while
the mesh, basis, weak-form, sparse assembly, and linear solve workflow remain
the same.

For example, the following assembles the Laplace operator for a unit square
with a uniform permittivity.  The left and right Dirichlet values can then be
handled by the application's preferred sparse-system utility.

```python
import numpy as np
import skfemntv
from skfemntv.helpers import ddot, grad

mesh = skfemntv.MeshTri.init_tensor(
    np.linspace(0., 1., 21), np.linspace(0., 1., 21)
)
basis = skfemntv.Basis(mesh, skfemntv.ElementTriP1())

@skfemntv.BilinearForm
def electrostatic_laplace(u, v, w):
    epsilon = 1.0
    return epsilon * ddot(grad(u), grad(v))

stiffness = skfemntv.asm(electrostatic_laplace, basis)
left = basis.get_dofs(
    nodes=lambda x: np.isclose(x[0], 0.)
).all()
right = basis.get_dofs(
    nodes=lambda x: np.isclose(x[0], 1.)
).all()
```

This is a compact example of extending the package beyond structural
mechanics: the physical meaning changes from displacement and stress to
potential and electric field `E = -grad(phi)`, but the readable Python weak
form stays intact.

The complete script is [`examples/plate_with_hole_tension.py`](examples/plate_with_hole_tension.py).
It uses PyVista for rendering.  Install the optional visualization packages
with `python -m pip install pyvista vtk`, then run
`python examples/plate_with_hole_tension.py` to regenerate the image.  The
mesh uses a polygonal approximation of the circular hole to keep mesh
generation dependency-free.  The script also reports the free-DOF residual
normalized by the applied boundary contribution, which is about `3e-14` for
this run.

For the supported mesh, element, basis, and assembly subset, the compatibility
target is stronger: changing

```python
import skfem
```

to

```python
import skfemntv as skfem
```

should preserve public DOF ordering and numerical results.  Higher-order DOFs
are numbered by topological entity, following scikit-fem: vertices first,
shared edges or facets next, and element interiors last.  This makes assembled
vectors and matrices directly interchangeable without a coordinate
permutation.

Solver conveniences such as `solve`, `enforce`, and `penalize` are outside the
native assembly backend.  Applications may continue to use scikit-fem or SciPy
for these operations.

## Capabilities

- functional, linear, bilinear, and cross-bilinear assembly
- meshes, elements, bases, tabulation, and quadrature
- caller-supplied tensor coefficient assembly
- threaded native kernels and sparse scatter
- cut-cell and implicit-interface quadrature
- interface supermeshes and contact-facet search

Scalar anisotropic diffusion uses the same readable helper syntax as
scikit-fem:

```python
from skfemntv.helpers import dot, grad, mul

@skfemntv.BilinearForm
def anisotropic_diffusion(u, v, w):
    return dot(mul(w.diffusion, grad(u)), grad(v))

A = np.array([[2.0, 0.4], [0.1, 0.8]])
matrix = skfemntv.asm(anisotropic_diffusion, basis, diffusion=A)
```

The recommended public layout is the scikit-fem component-first form
`(dim, dim, entities, quadrature)`.  A constant `(dim, dim)` tensor is also
accepted and broadcast over all integration points.  The entity-first native
layout `(entities, quadrature, dim, dim)` remains accepted for low-level
compatibility, but is not the preferred user-facing representation.  The
native kernel uses that internal layout so each quadrature-local tensor is
contiguous in memory.  Nonsymmetric tensors are retained rather than silently
symmetrized.  This first typed operation targets scalar H1 fields;
vector-field fourth-order constitutive tensors remain separate work.

For numeric component-first arrays, `skfemntv.helpers` also provides the
scikit-fem-compatible `det`, `inv`, `eye`, and `cross` operations. They work on
constant tensors and quadrature data. Using these operations on symbolic
nonlinear form expressions is currently rejected by the typed form compiler.

Coefficient components follow scikit-fem's component-first convention as
well: `w.material[0]` selects axis 0.  Multiple named fields can be used in a
single typed form, for example:

```python
@skfemntv.BilinearForm
def material_form(u, v, w):
    return (
        w.material[0] * dot(u, v)
        + w.material[1] * ddot(grad(u), grad(v))
        + dot(mul(w.diffusion, grad(u)), grad(v))
    )
```

Component indices must be integers.  Missing fields and out-of-range
components are reported with the coefficient name before assembly.

The mathematical compatibility boundary is also available at runtime:

```python
skfemntv.supports("space.h1")       # True
skfemntv.supports("space.hcurl")    # False (experimental opt-in required)
skfemntv.require_capability("space.hcurl")  # precise capability error
```

Use `skfemntv.capabilities()` for the complete machine-readable registry.
Experimental capabilities require `include_experimental=True`; solver policy
is explicitly marked external.

`Basis.project` provides experimental L2 projection for scalar and vector
P0/P1/P2 bases on affine triangular and tetrahedral meshes. It currently
projects over the full basis; composite fields and curved meshes are not
supported.

### Experimental H(curl) simplex slices

The lowest-order first-family Nédélec element is available as a deliberately
small experimental API for affine triangular meshes:

```python
import numpy as np
from scipy.sparse.linalg import spsolve
import skfemntv

mesh = skfemntv.MeshTri.init_tensor(
    np.linspace(0.0, 1.0, 9), np.linspace(0.0, 1.0, 7)
)
basis = skfemntv.AffineTriN1Basis(mesh, intorder=3)
assembler = skfemntv.TriN1Assembler(basis)
linear_assembler = skfemntv.TriN1LinearAssembler(basis)
matrix = assembler.assemble_maxwell(
    mass_coefficient=1.0, curl_coefficient=0.05
).copy()

boundary = basis.boundary_dofs()
free = np.setdiff1d(np.arange(basis.N), boundary)
load = linear_assembler.assemble_vector_load(
    lambda x: np.array((0.0 * x[0], np.sin(np.pi * x[0])))
).copy()
solution = np.zeros(basis.N)
solution[free] = spsolve(matrix[free][:, free], load[free])
```

Analytic vector fields can be interpolated by their globally oriented edge
moments and evaluated at basis quadrature points:

```python
def field(x):
    return np.array((0.0 * x[0], x[0] * x[1]))

edge_values = basis.interpolate_edge_moments(field)
values = basis.evaluate(edge_values)       # (component, cell, quadrature)
curls = basis.evaluate_curl(edge_values)   # (cell, quadrature)
points = basis.global_coordinates          # (component, cell, quadrature)
```

The tetrahedral counterpart uses the same explicit workflow:

```python
axis = np.linspace(0.0, 1.0, 5)
mesh = skfemntv.MeshTet.init_tensor(axis, axis, axis)
basis = skfemntv.AffineTetN1Basis(mesh, intorder=3)
assembler = skfemntv.TetN1Assembler(basis)
linear_assembler = skfemntv.TetN1LinearAssembler(basis)

# TetN1 integration is cell-parallel; applications retain thread ownership.
with skfemntv.thread_limit(4):
    matrix = assembler.assemble_maxwell(
        mass_coefficient=1.0, curl_coefficient=0.05
    ).copy()

load = linear_assembler.assemble_vector_load(
    lambda x: np.array((1.0 + 0.0 * x[0], x[0], -x[1]))
).copy()
boundary = basis.boundary_dofs()
```

TetN1 values and vector curls both use
`(basis, component, cell, quadrature)` publicly; evaluated fields and curls
use `(component, cell, quadrature)`.  The complete runnable example is
`examples/hcurl_tet_n1_maxwell.py`.

These slices include oriented global edge numbering, affine covariant Piola
mapping, mass/curl-curl/Maxwell CSR assembly, coefficient fields, memory
preflight, and boundary-edge selection.  They support only affine triangles and
tetrahedra with one lowest-order tangential moment per edge.  Neither basis is
accepted by the general `Basis` or `asm`; curved mappings, higher order, and
solver policy are not provided.  Both bases support edge-moment interpolation
and quadrature-point evaluation.  The triangle example is
`examples/hcurl_tri_n1_maxwell.py`.

`TriN1Assembler` uses NumPy/SciPy integration.  `TetN1Assembler` uses a fused
cell-parallel C++ integration kernel and the same fixed-CSR ownership model.
Names describe their mathematical spaces rather than implementation details.
Assembly methods reuse
and overwrite one result object, so call `.copy()` before another assembly when
results must be retained.  Public basis values use component-first shape
`(basis, component, cell, quadrature)` and curls use
`(basis, cell, quadrature)` for triangles and
`(basis, component, cell, quadrature)` for tetrahedra; coefficient fields use
`(cell, quadrature)`.

`basis.geometry_diagnostics` reports minimum signed/absolute Jacobian
determinants, minimum cell area/volume, an edge-based aspect indicator, and the
number of negatively oriented cells.  An optional `max_aspect_ratio` constructor
argument rejects meshes beyond an application-selected quality limit.  Tests
cover distorted meshes and mixed cell orientations, including manufactured
PDE convergence.  With this sharply limited scope, `space.hcurl` is declared
experimental; this does not imply support for general meshes or form syntax.

### Assembly memory preflight

Estimate native bilinear-assembler memory before its CSR pattern and scatter
map are allocated:

```python
estimate = skfemntv.estimate_bilinear_memory(basis)

print(skfemntv.format_bytes(
    estimate.construction_peak_total_bytes_upper_bound
))

if not estimate.fits_in(available_bytes, safety_factor=1.25):
    raise MemoryError("assembly does not fit the configured memory budget")
```

The same check can be enforced atomically before native CSR and scatter-map
allocation:

```python
assembler = skfemntv.NativeBilinearForm(
    basis,
    memory_limit_bytes=16 * 1024**3,
    memory_safety_factor=1.25,
)
```

`NativeCrossBilinearForm` accepts the same options.  A failed check raises
`AssemblyMemoryBudgetError` containing the estimate, configured budget,
safety-adjusted requirement, and largest estimated allocation.  The default
safety factor is 1.25.  Supplying no limit preserves existing behavior and does
not impose a process-wide policy.  Segmented `CutCellBasis` bilinear and cross
assemblers use active-cell counts and flattened cut quadrature in the same
preflight and budget API.

The report separates memory already retained by the basis from persistent
native assembler allocations and temporary CSR-pattern construction memory.
Its CSR nonzero count is a conservative, allocation-free upper bound.  Python
object overhead, allocator fragmentation, and coefficient arrays supplied at
assembly time are intentionally excluded and listed in the report assumptions.

Use `estimate_bilinear_memory(test_basis, trial_basis)` for a rectangular
cross-bilinear block.  For mixed spaces,
`estimate_composite_bilinear_memory(basis, field_pairs=...)` estimates the
selected independently cached blocks; omitting `field_pairs` assumes all
ordered field pairs are assembled once.  Repeat a field pair when distinct
value/gradient contraction kinds create multiple cached native blocks.

## Assembly performance

The following reference run compares warm-cache Poisson P1 matrix and
right-hand-side assembly using identical meshes and spaces.  Results depend on
hardware, thread affinity, and package versions; they are an example rather
than a universal performance guarantee.

![Poisson assembly scaling](https://raw.githubusercontent.com/kevin-tofu/scikit-fem-native/main/benchmarks/compare-with-skfem/results/poisson-linux-x86_64.png)

See the [benchmark methodology](benchmarks/compare-with-skfem/README.md) and
[recorded environment and values](benchmarks/compare-with-skfem/results/poisson-linux-x86_64.md)
for details.

See [DEVELOPMENT.md](DEVELOPMENT.md) for the motivation, architecture boundary,
and development principles.

## Documentation

Read the [scikit-fem-native documentation](https://scikit-fem-native.readthedocs.io/en/latest/),
including the precise
[scikit-fem compatibility boundary](https://scikit-fem-native.readthedocs.io/en/latest/scikit-fem-compatibility.html)
and the
[assembly performance methodology](https://scikit-fem-native.readthedocs.io/en/latest/scikit-fem-performance.html).

Build the Sphinx documentation locally:

```bash
python -m pip install -e '.[docs]'
sphinx-build -M html docs docs/_build
```

Open `docs/_build/html/index.html` after the build completes.

## Development

```bash
python -m pip install -e '.[test]'
pytest -q
```
