geossm.covmodel.buildMesh2d_density

geossm.covmodel.buildMesh2d_density(points, lowrank, domain=None, max_edge=None, min_edge=None, offset=None, cutoff=None, min_angle=21.0, snap_to_points=True, relax_iters=2, tol=0.05, max_iter=20, seed=None)

Build a 2D triangular mesh whose interior vertices are a density- matched subset of points, rather than free vertices placed by a smooth size field (as buildMesh2d does).

Where buildMesh2d steers a Delaunay mesher’s local element size to approximate the density of points, this function instead: (1) picks landmark locations via k-means on points – since k-means codebook density tracks source density, this directly matches mesh-vertex density to point density; (2) by default snaps each landmark to its nearest not-yet-used point (a greedy closest-pair matching), so landmarks coincide exactly with observed locations instead of merely being near them, letting the mesh capture the observed variability directly at those nodes; (3) embeds the landmarks into the gmsh mesh as fixed points – gmsh’s optimizer moves every other node but never relocates or removes an embedded point, so the overlap from (2) survives mesh generation and smoothing exactly. Between landmarks, gmsh fills in the rest from their own local spacing (Mesh.MeshSizeFromPoints), which lets element size track landmark density rather than being held to a single max_edge everywhere.

That said, max_edge is still enforced as a hard coverage bound, not just a size cap where landmarks are dense: any location more than max_edge from every landmark – the outer offset buffer (which has no landmarks at all), a bay or strait inside the convex hull but outside a concave or multi-part domain (e.g. a country plus its islands), or, just as much, a genuinely sparse or snap_to_points- relaxation-drifted pocket inside the interest domain – gets an extra lattice (“filler”) vertex there too, embedded exactly like a landmark. This matters because an unbounded gap doesn’t just mean a coarse patch of mesh: it means an unbounded-variance patch of the SPDE/FEM latent field, which a data-density-matched mesh alone cannot prevent if points happens to be sparse near, say, the domain’s edge – a common real-world pattern. (A runtime “outside the interest domain -> cap at max_edge” size field, and later a version that only filled outside the interest domain, were both tried before this lattice; the field’s discontinuity made gmsh’s mesher visibly struggle on a finely-detailed real coastline, and restricting filler to outside the interest domain left real gaps wherever points was itself sparse.) The lattice is pulled back about half a cell from the outer boundary, rather than filling all the way up to it: a raw grid clipped flush against a slanted or jagged edge produces its own near- zero-angle slivers, so that margin is instead left to gmsh’s ordinary boundary-conforming triangulation, which bridges it cleanly from the boundary vertices already there. Either way, max_edge (via the filler lattice) and the landmarks’ own local spacing are what keep interior angles away from zero: an abrupt jump from a tight cluster of landmarks straight to a coarse neighboring region would force sliver triangles, so gmsh grades a few extra vertices in between as needed.

Because that grading – now including the coverage-driven filler – adds vertices beyond the landmarks themselves, asking k-means for exactly round(lowrank * len(points)) landmarks can overshoot the vertex-count target, sometimes substantially if lowrank is small relative to how large the interest domain is compared to max_edge (the filler coverage floor doesn’t shrink just because lowrank asked for fewer landmarks). So, mirroring buildMesh2d’s bisection over a size-field scale, this function bisects over the number of landmarks requested instead, regenerating the mesh each time, until the total vertex count inside the interest domain (landmarks plus grading) is within tol of the target – see buildMesh2d’s lowrank for what “interest domain” means here. When that floor exceeds the target no amount of bisecting can close the gap, and a warning says so; a coarser max_edge is the way to actually lower the floor, not a smaller lowrank.

Landmarks closer together than cutoff are merged, which trades off against min_angle: forcing two near-duplicate data points to both remain exact, fixed vertices is often what produces a sliver triangle that no amount of optimization can fix, since fixed points cannot be moved. Widen cutoff (or reduce lowrank) if the min_angle warning fires often.

Parameters:

• points ((n, 2) array_like) – Observed locations.

• lowrank (float) – Value in [0, 1]. Target number of vertices inside the interest domain, as a fraction of len(points) – see above. Values that would ask for fewer than 3 landmarks, including 0, are floored to 3 (a single triangle), the smallest possible triangulation.

• domain (Polygon, MultiPolygon, or (possibly nested) list of these, optional) – The scientific-interest domain; see buildMesh2d. Defaults to the convex hull of points.

• max_edge (float, optional) – Largest allowed triangle edge length: a coverage bound enforced by the filler lattice anywhere a location is more than max_edge from every landmark, inside the interest domain as much as outside it, and an upper bound on each landmark’s own local- spacing size – see above. This can conflict with a small lowrank on a large interest domain, since it sets a vertex-count floor lowrank cannot go below; a warning says so rather than silently leaving the mesh under-covered. Defaults to 1/15 of the domain’s bounding-box diagonal.

• min_edge (float, optional) – Smallest allowed triangle edge length. Defaults to max_edge / 10.

• offset (float, optional) – Buffer added around domain so the mesh extends past the data. Defaults to max_edge.

• cutoff (float, optional) – Minimum allowed separation between landmarks, and independently between a landmark and any filler vertex (see above) – a mesh- quality knob, unlike buildMesh2d’s cutoff: it never merges or discards any of points itself (only exact duplicates are merged, always, regardless of this setting), so it can never make lowrank (e.g. lowrank=1, “use every point”) see fewer points than you passed in. Left at its default, this is derived from target_n and the interest domain’s area instead of max_edge (unlike buildMesh2d): a large max_edge, appropriate for a sparse, data-free outer buffer, has nothing to do with how tightly target_n landmarks need to pack, and a max_edge/5 default would otherwise silently cap the achievable landmark count well below target_n. That default is also shrunk automatically (up to a few times) if it still leaves too few landmarks to reach target_n. Passing cutoff explicitly disables both and is honored exactly for landmark/filler spacing, even if that undershoots target_n (reported via a warning).

• min_angle (float, default 21.0) – Target minimum interior angle (degrees). This is a soft target graded around the fixed landmarks; it is not enforced on the landmarks’ own placement, which is data-driven, not quality-driven – see above.

• snap_to_points (bool, default True) – If True (recommended), landmarks coincide exactly with observed points (greedy nearest-pair matching from k-means centroids to points). If False, landmarks are left at the k-means centroids (relaxed by relax_iters) – close to the data but generally not exactly on it, which can give a smoother, better-quality triangulation.

• relax_iters (int, default 2) – Only used when snap_to_points=False. Number of Delaunay-neighbor averaging passes applied to the k-means centroids before they are embedded, nudging landmarks that ended up too close to a neighbor toward more even local spacing. Centroids alone are not enough: when target_n (via lowrank) approaches len(points), k-means with that many clusters degenerates to one point per cluster, so every centroid lands exactly on its point regardless of snap_to_points – this is what actually gives snap_to_points room to move landmarks off the raw data in that regime. Set to 0 to disable and use the raw (possibly point-coincident) centroids.

• tol (float, default 0.05) – Relative tolerance on the vertex-count target used to stop the landmark-count search. Landmark placement (via k-means) and mesh grading are both a little noisy as the requested count changes by one, so this defaults looser than buildMesh2d’s tol.

• max_iter (int, default 20) – Maximum number of mesh (re)generations used by the search.

• seed (int, optional) – Seed for the k-means initialization. Also fixes it internally across the search’s repeated k-means calls (even if left None, in which case a seed is drawn once and reused for this call only), so that requesting fewer landmarks is what changes the mesh, not fresh k-means randomness.

Returns:

• mesh (meshio.Mesh)

• convex_hull (shapely.geometry.Polygon) – The (buffered) convex hull the mesh was built over – domain (or its default) widened to its convex hull and extended by offset.