From RVE to Mesh: A Pipeline for Heterogeneous Continua

Two earlier posts on this site developed complementary halves of a single picture. Reservoir Geometry took an anisotropic permeability field as given and built a Riemannian manifold on top of it with . Numerical Analysis via Discrete Exterior Calculus took a manifold as given and discretised its Laplace-Beltrami operator using the Hodge star on a simplicial mesh. This post closes the loop on both ends. Upstream of the reservoir post, we ask where comes from, and answer with mean-field homogenisation on a representative volume element. Downstream of the DEC post, we ask how the mesh should adapt to the resulting geometry, and answer with curvature-driven remeshing using the closed-form scalar curvature .

The result is a single arc from grain-scale microstructure to a converged finite-dimensional solution: microstructure feeds homogenisation, homogenisation produces a tensor field, the tensor field is a metric, the metric induces a Hodge star, the Hodge star discretises the Laplacian, and the curvature of the metric tells the mesh where to refine. Every arrow in that chain is a mathematical operation with a closed-form or algorithmic expression, and every arrow is implemented in the cartan workspace.

1. Scale separation and the representative volume element

A porous or fractured rock is heterogeneous on the scale of pores and grains, with length scale on the order of micrometres, and it is homogeneous on the scale of reservoirs, with length scale on the order of hundreds of metres. A representative volume element (RVE) is a sub-volume of the rock at an intermediate scale , satisfying

The RVE is large enough that the statistics of its microstructure (pore sizes, fracture orientations, grain packing) match those of the bulk rock, and small enough that macroscopic fields (pressure, stress) are effectively constant across it. The existence of such an intermediate scale is the scale separation hypothesis. It is the assumption that licenses every homogenisation scheme discussed below.

Introduce the small parameter and the fast variable . A heterogeneous tensor field oscillating on the RVE scale has the form , where is -periodic in . The goal of homogenisation is to produce a slowly varying effective tensor that captures the macroscopic behaviour of as .

2. Two-scale asymptotic homogenisation

The rigorous route to is two-scale asymptotic expansion, developed for periodic structures by Bensoussan, Lions and Papanicolaou [5] and placed on the footing of two-scale convergence by Allaire [6]. Consider the Darcy pressure equation on the reservoir, and expand the solution in powers of :

with each assumed -periodic in . The spatial derivative splits according to the chain rule: , acting on functions of both slow and fast variables. Substituting into the Darcy equation and collecting powers of produces a hierarchy of problems.

The balance forces to be independent of the fast variable, . The balance gives the cell problem: find , periodic in , such that

The cell problem is a linear elliptic boundary-value problem on the periodic unit cell , parameterised by the slow variable and solved once per macroscopic point (or once per RVE if the microstructure is statistically homogeneous). Its solution is the corrector that encodes the fine-scale response to a unit macroscopic gradient in direction .

The balance yields the homogenised macroscopic equation , with

The cell average is taken with respect to the fast variable; the slow variable is a parameter. Symmetry and positive definiteness of follow from symmetry and positive definiteness of the microstructural tensor , so at every .

3. Mean-field homogenisation and the Eshelby localisation

Two-scale asymptotics requires solving the cell problem numerically on each RVE. For fractured and porous media with a small number of inclusion families (matrix, pores, cracks with a few orientation classes), a closed-form alternative is available: mean-field homogenisation in the Zaoui [2] and Ponte Castañeda-Willis [4] tradition. The starting point is the Eshelby inclusion problem [3].

For a sphere in an isotropic matrix (, ), the Hill tensor is . Closed forms for general spheroids in isotropic and transversely isotropic matrices are tabulated in the ECHOES documentation [10] and in Milton's monograph [1].

From Eshelby's theorem one constructs, for each inclusion family with volume fraction and conductivity contrast , a dilute localisation tensor

The effective conductivity is then a weighted combination of the . Three schemes cover most use cases.

Dilute scheme. Non-interacting inclusions in a host matrix:

Mori-Tanaka scheme. Each inclusion sees the matrix at the average field of the composite, enforcing that concentration tensors sum to the identity:

Self-consistent scheme. Each phase is embedded in the effective medium itself, yielding an implicit equation solved by fixed-point iteration:

Each scheme produces an SPD tensor field at every macroscopic point. For a fractured reservoir whose crack statistics (density, orientation distribution, aspect ratio) vary with , Mori-Tanaka gives the effective permeability field the reservoir post takes as input.

4. The metric and its fibre bundle

Given at every , define the Riemannian metric

This is the construction developed in the reservoir geometry post. The reservoir becomes a Riemannian manifold , streamlines are geodesics, pressure diffusion is governed by the Laplace-Beltrami operator, and permeability contrasts are refractive discontinuities.

The structure is cleaner when viewed as a fibre bundle. Let be the bundle whose fibre at is , the space of symmetric positive definite matrices. The effective permeability is a section of this bundle,

Two operations sit naturally on this bundle.

Base-horizontal operations are calculus along with the metric at each point: the Levi-Civita connection, the Laplace-Beltrami operator, geodesic streamlines, scalar curvature. These are the content of the reservoir post.

Fibre-vertical operations are geometry inside at a fixed base point: the affine-invariant Riemannian metric on [9], geodesic interpolation between measured values at neighbouring wells, the Fréchet mean of a distribution of SPD samples. These are used when is known only at discrete measurement points and must be interpolated to the rest of .

The same structure that parameterises the fibre is used by every stage of the pipeline: homogenisation outputs a point in per RVE, the metric inverts that point at each , the Hodge star (next section) is parameterised by the same tensor, and mesh quality at each triangle is measured against the same inner product.

5. The Hodge star from the homogenised metric

The DEC post developed discrete exterior calculus on a flat Euclidean background, following the framework of Hirani [7] and Desbrun, Hirani, Leok and Marsden [8]. The extension to a Riemannian manifold is where the metric enters the discretisation. In the continuous setting, the Hodge star is defined by

where the pointwise inner product on -forms is constructed from (and its inverse ), and is the Riemannian volume form. The Laplace-Beltrami operator on 0-forms is .

On a simplicial mesh of , the Hodge star becomes a diagonal (or block-diagonal) linear map between cochain spaces. For the scalar Poisson problem on vertices, the diagonal circumcentric Hodge star maps a primal edge of length (measured in the metric ) to its dual edge of length , with ratio

Every geometric quantity in [eq:hodge-diagonal] is computed in the metric induced by . The length of a primal edge from vertex to vertex is

with a geodesic in . For a fine mesh, the straight-line approximation

with the midpoint is second-order accurate in the mesh size, which matches the order of the linear finite-element discretisation used downstream.

The Laplace-Beltrami operator assembled from is the discrete analogue of the continuous operator from the reservoir post. Its null space consists of constants (for Neumann boundary conditions), its spectrum depends on , and its condition number is governed by the anisotropy ratio of . Every property of that matters for pressure diffusion is inherited through the Hodge star.

6. Scalar curvature as an a priori refinement indicator

The closed-form scalar curvature of the isotropic permeability metric is the capstone identity of this pipeline: it is the single formula that starts from the homogenisation output and ends at the refinement indicator used by cartan-remesh. Every earlier section feeds into it, and every later step consumes it.

The formula [eq:R-scalar] is the hinge of the pipeline. Its inputs are produced by homogenisation (section 3); its output drives the refinement criterion below. Two scalar derivatives of the homogenised permeability are everything the remesher needs.

Regions of high are regions where the metric is varying rapidly, which are exactly the regions where streamlines bend sharply, Jacobi fields spread or focus, and pressure transients accumulate error. These are the regions where a finite-element mesh must be refined if the discrete solution is to converge uniformly.

Classical adaptive finite-element codes use a posteriori residual estimators: solve on a coarse mesh, estimate the local error from the residual, refine, repeat. The approach works but requires a full solve per iteration. A complementary strategy is available once the metric is known analytically: use itself as an a priori refinement indicator, evaluated on vertices of the current mesh before any solve.

The criterion [eq:curvature-cfl-criterion] is the geometric analogue of a CFL condition: it balances the local edge length against the local curvature scale, so that no edge spans more than units of "curvature-length." The cartan-remesh crate implements exactly this driver: adaptive_remesh applies split, collapse, flip, and shift passes based on a RemeshConfig whose curvature field is the pointwise .

For anisotropic , equation [eq:R-scalar] generalises. The scalar curvature of on is a contraction of second derivatives of with , computable symbolically once and for all. For the transversely isotropic reservoir model , the expression is tractable in closed form; for fully anisotropic , it is evaluated numerically from and its first two derivatives using finite differences on the mesh. Either way, the cost per vertex is , versus the cost of an a posteriori residual estimate.

7. The full pipeline

Assembling the pieces gives a seven-stage pipeline from microstructure to discrete solution.

Each stage is a named module in the cartan workspace. Stage 2 is scheduled for cartan-core's homogenisation submodule; stage 3 is trivially the inverse in SPD(3); stages 4-6 are already in cartan-dec and cartan-remesh; stage 7 uses cartan-dec's sparse solvers (or an external solver via the Rust bindings).

8. Worked example: a transversely isotropic fractured reservoir

As a minimal end-to-end instance, consider a fractured sandstone with matrix permeability and a single family of horizontal penny-shaped cracks with aspect ratio , density , and orientation normal to . The crack density varies with depth on a characteristic scale .

Stage 1. The microstructural data is .

Stage 2. For penny-shaped cracks, the Mori-Tanaka effective permeability is transversely isotropic, , with

in the dilute limit (corrections of order from full Mori-Tanaka). The horizontal permeability oscillates with depth; the vertical permeability is unchanged because cracks normal to conduct horizontally only.

Stage 3. The metric is , varying with depth.

Stage 5. The Laplace-Beltrami operator on vertex 0-forms is

(with the standard chain-rule factors). The discrete form uses diagonal Hodge stars whose entries are or depending on edge orientation.

Stage 6. The scalar curvature picks up terms. Where varies fastest (near the zero-crossings of ), is largest, and the curvature-CFL criterion refines the mesh in horizontal slabs tracking the inflection surfaces. The resulting mesh has high density in the transition layers and low density in the plateaus, matching the -optimal node distribution for the pressure problem up to a constant.

The pipeline produces, from four numbers and one shape parameter , a converged discrete pressure solution whose mesh reflects the depth structure of the fracture field. No empirical error estimator, no residual-based adaptivity, no user-tuned refinement zones.

9. Discussion

Why the pipeline is closed. Each stage consumes the output of the previous one through a mathematical operation, not an engineering convention. Homogenisation produces a tensor of a prescribed algebraic type (symmetric positive definite). That type is exactly what a Riemannian metric needs. The metric is exactly what a Hodge star needs. The Hodge star is exactly what a Laplace-Beltrami operator needs. The Laplacian's coefficients are exactly what a scalar curvature formula needs. The curvature is exactly what a mesh adaptation driver needs. No stage invents structure that the next stage ignores, and no stage asks for structure that the previous stage failed to provide.

Why threads through everything. Symmetric positive definite matrices are the right algebraic object for linear anisotropic response in any setting where energy is positive. The space is a Riemannian manifold in its own right, with the affine-invariant metric that respects congruence. When this structure is respected (geodesic interpolation in the fibre, affine-invariant distances for mesh comparison, Fréchet means for uncertainty quantification), every downstream operation inherits good properties: positive-definiteness preservation, invariance under coordinate rotation, stable interpolation under rank-one updates. Violating the SPD structure (Euclidean interpolation of entries, for example) is the standard way that adaptive reservoir codes produce non-physical negative eigenvalues under mesh refinement.

What is not in this post. Nonlinear homogenisation (where depends on the local field), time-dependent homogenisation (ageing viscoelasticity, where the Hill tensor acquires a memory kernel), and stochastic homogenisation (random microstructure, where is itself a random tensor with a limiting distribution) all sit outside the static linear pipeline described here. Each extension preserves the backbone of the present picture and adds one layer: a fixed-point iteration for nonlinear, a convolution in time for viscoelastic, a mean and covariance on for stochastic. The cartan workspace is structured to accommodate these as future crates without disturbing the flat linear case.

Connection to the wider program. The RVE-to-mesh pipeline is one instance of a pattern that recurs in reservoir geometry, DEC numerical analysis, and the active-matter and membrane work developed elsewhere on this site. In each case a heterogeneous physical system is described by a tensor field on a base manifold; the tensor field is interpreted as a geometric object (metric, connection, curvature); that geometry is discretised by DEC; and the discretisation is refined adaptively using an analytic curvature quantity. The reservoir is one case study. The machinery generalises.

References

1. G. W. Milton. The Theory of Composites. Cambridge University Press, 2002.
2. A. Zaoui. "Continuum micromechanics: survey." Journal of Engineering Mechanics, 128(8):808–816, 2002.
3. J. D. Eshelby. "The determination of the elastic field of an ellipsoidal inclusion, and related problems." Proceedings of the Royal Society A, 241:376–396, 1957.
4. P. Ponte Castañeda and J. R. Willis. "The effect of spatial distribution on the effective behavior of composite materials and cracked media." Journal of the Mechanics and Physics of Solids, 43(12):1919–1951, 1995.
5. A. Bensoussan, J.-L. Lions, and G. Papanicolaou. Asymptotic Analysis for Periodic Structures. North-Holland, Amsterdam, 1978.
6. G. Allaire. "Homogenization and two-scale convergence." SIAM Journal on Mathematical Analysis, 23(6):1482–1518, 1992.
7. A. N. Hirani. Discrete Exterior Calculus. PhD thesis, California Institute of Technology, 2003.
8. M. Desbrun, A. N. Hirani, M. Leok, and J. E. Marsden. "Discrete exterior calculus." arXiv:math/0508341, 2005.
9. X. Pennec, P. Fillard, and N. Ayache. "A Riemannian framework for tensor computing." International Journal of Computer Vision, 66(1):41–66, 2006.
10. J.-F. Barthélémy. ECHOES: Extended Calculator of HOmogEnization Schemes. jfbarthelemy.github.io/echoes, 2024.