geossm.covmodel.FEMSolver

class geossm.covmodel.FEMSolver(meshio_obj, domain=None, verbose=True, stats=True)[source]

Bases: object

Low-level finite-element engine behind spdeAppoxCov.

Wraps a triangular mesh (e.g. built by buildMesh2d) in an MFEM order-1 \(H^1\) finite element space, and assembles the sparse mass and stiffness matrices that spdeAppoxCov.precision combines into the SPDE precision matrix. It also classifies mesh vertices as “inner” (covered by domain, the region of scientific interest) or “outer” (part of the surrounding buffer, whose role is only to push the FEM’s Neumann boundary effects away from domain), evaluates basis functions at arbitrary points (getBasis), and reports mesh quality diagnostics (compute_stats, summary).

Most users go through spdeAppoxCov instead, which builds and owns a FEMSolver internally in spdeAppoxCov.setup(). FEMSolver is exposed directly (e.g. via spdeAppoxCov.fem_solver) for lower-level access to the mesh, the assembled matrices, and mesh diagnostics.

Notes

Internally, the input meshio.Mesh is converted to an mfem.Mesh via meshio_to_mfem_mesh (2D triangular meshes only); vertices and elements are kept in the order given by the mesh, with no reordering or optimisation performed by this class.

Methods

compute_angles

The three interior angles (degrees) of a single triangle.

compute_area

Area of a single triangle (shoelace formula).

compute_stats

Interior angles and areas of every mesh triangle (via compute_angles/compute_area); used by __init__ (when stats=True) to populate angles/areas and warn about poor-quality (small-angle) meshes.

covers_domain

Check that the mesh covers the FEM computational domain Ω.

distance

Alias for get_distance.

generate_summary

Build the (label, value) rows shown by summary, as two lists of (str, [str]) pairs (left and right columns).

generate_summary_header

Cheap mesh-size rows (nelements/nbElements, direct O(1) MFEM accessor calls) - unlike generate_summary(), does not touch nvertex/n_inner_points/n_outer_points/angles/areas, which involve a Python loop over every mesh element or a reduction over the (possibly large) angle/area arrays.

getBasis

Evaluate the FE basis functions at a set of physical points.

get_distance

Flat Euclidean pairwise distance between points and the mesh vertex coordinates (or between mesh vertices and themselves).

isinner

Classify vertices as inner or outer with respect to domain.

plot_mesh

Plot the mesh: all triangle edges, inner/outer vertices, boundary edges (dashed red), and the domain polygon(s) (orange outline).

print_info

Print msg prefixed with a UTC timestamp.

summary

Return a statsmodels.iolib.summary.Summary table of mesh/FE diagnostics (vertex/element counts, angle & area ranges, DOFs, matrix shapes, ...); also backs __str__.

Attributes

GetVSize

Total number of vector DOFs in fespace (ndofs * fespace_order).

angles

(nelements, 3) ndarray of the triangles' interior angles (degrees), or None if __init__ was called with stats=False.

areas

(nelements,) ndarray of triangle areas, or None if __init__ was called with stats=False.

boundary_vertex

(nbElements, 2) ndarray of int -- vertex-index pairs of the boundary elements (see nbElements).

box

Mesh bounding box as [minx, miny, maxx, maxy] (rounded to 2 decimals).

domain

region(s) of scientific interest used to classify vertices as inner/outer (see _validate_domain / isinner).

domain_hull

Convex hull of the union of domain -- see _domain_hull.

effective_dofs

Return the effective number of DOFs (associated with vertices)

elements

(nelements, 3) ndarray of int -- vertex indices (into vertex) making up each mesh triangle.

fespace

Return the finite elment space

fespace_order

Vector dimension of fespace (1 for the scalar H1 space used here).

geoshape

MFEM geometry type of the mesh elements (e.g. triangle).

getBoundaryEdge

(nbElements, 2) ndarray of int -- vertex-index pairs for each boundary edge.

inner

Boolean array, one entry per vertex, True where the vertex is classified as inner (covered by domain); see isinner.

inner_points

Return the point of the grid

mass

Sparse (diagonal, lumped) mass matrix \(C\), shape (effective_dofs, effective_dofs).

mesh

The underlying mfem.Mesh (converted from the input meshio.Mesh).

n_inner_points

Return the point of the grid

n_outer_points

Return the number of the boundary point of the grid

n_totpoints

Return the total point number of the grid

nbEdges

Total number of mesh edges.

nbElements

Number of boundary elements (edges, for a 2D mesh).

ndofs

Number of scalar degrees of freedom in fespace (all mesh vertices, before cutting to effective_dofs).

nelements

Get the number of discretisation elements (like the number of triangle)

nvertex

Number of mesh vertices (= effective_dofs).

outer_points

Return the boundary point of the grid

shape

Tuple (nvertex, nelements, nbElements).

stiff

Sparse stiffness matrix \(G\) (the <grad(phi_i), grad(phi_j)> FEM Laplacian term), shape (effective_dofs, effective_dofs).

totpoints

Return the total point (inner + boundary) of the grid (same as vertex)

vertex

Returns the coordinates of mesh vertices that have a DOF associated with them in the given H1_FECollection(1, dim) finite element space.

property mesh

The underlying mfem.Mesh (converted from the input meshio.Mesh).

property inner

Boolean array, one entry per vertex, True where the vertex is classified as inner (covered by domain); see isinner.

isinner(points=None)[source]

Classify vertices as inner or outer with respect to domain.

A point is inner if it is covered by any single polygon in domain – each polygon is tested on its own (not on the convex hull of their union), so for a multi-polygon domain (e.g. Sicily, Sardinia, mainland Italy) a mesh vertex sitting between two landmasses, inside their shared convex hull but outside every individual polygon, is classified as outer.

getBasis(phy_points=None)[source]

Evaluate the FE basis functions at a set of physical points.

For each point, locates the mesh element that contains it and evaluates every basis function of that element there, giving a sparse “observation” matrix H such that H @ field_at_vertices interpolates the FEM field at phy_points.

Parameters:

phy_points ((p, 2) array_like, optional) – Physical-space points to evaluate the basis at. Defaults to the mesh vertices (vertex), giving the identity-like mapping used e.g. to sanity-check the FE space.

Returns:

  • count (int) – Number of points for which a containing mesh element was found (see mfem.Mesh.FindPoints).

  • notfindInx (ndarray of int) – Indices (into phy_points) of points for which no containing element could be found – point search is not guaranteed to succeed even for points that do lie inside the mesh.

  • H (scipy.sparse.csr_matrix, shape (p, effective_dofs)) – Basis functions evaluated at each point, one row per point and one column per (inner + outer) mesh vertex/DOF. Rows for points in notfindInx are all zero. Entries below 1e-8 are zeroed out for numerical stability.

plot_mesh(ax=None, figsize=(10, 8), title='Title', alpha_vertex=1, alpha_triangle=0.5, alpha_border=0.5)[source]

Plot the mesh: all triangle edges, inner/outer vertices, boundary edges (dashed red), and the domain polygon(s) (orange outline).

Parameters:
  • ax (matplotlib.axes.Axes, optional) – Axes to draw on. A new figure/axes is created if not given.

  • figsize (tuple of float, default (10, 8)) – Figure size, only used when ax is not given.

  • title (str, default "Title") – Plot title.

  • alpha_vertex (float, default 1) – Transparency of the inner/outer vertex markers.

  • alpha_triangle (float, optional) – Unused (kept for interface compatibility); triangle edges are always drawn with a fixed alpha of 0.5.

  • alpha_border (float, default 0.5) – Transparency of the dashed boundary-edge lines.

Returns:

matplotlib.axes.Axes

get_distance(points=None)[source]

Flat Euclidean pairwise distance between points and the mesh vertex coordinates (or between mesh vertices and themselves).

Parameters:

points ((p, 2) array_like, optional) – Query points. Defaults to vertex, giving the full vertex-to-vertex distance matrix.

Returns:

ndarray, shape (nvertex, nvertex) or (p, nvertex)

distance(points=None)[source]

Alias for get_distance.

property vertex

Returns the coordinates of mesh vertices that have a DOF associated with them in the given H1_FECollection(1, dim) finite element space.

property nvertex

Number of mesh vertices (= effective_dofs).

property nelements

Get the number of discretisation elements (like the number of triangle)

property elements

(nelements, 3) ndarray of int – vertex indices (into vertex) making up each mesh triangle.

property box

Mesh bounding box as [minx, miny, maxx, maxy] (rounded to 2 decimals).

property domain

region(s) of scientific interest used to classify vertices as inner/outer (see _validate_domain / isinner).

Type:

Tuple of shapely.geometry.Polygon

property domain_hull

Convex hull of the union of domain – see _domain_hull.

covers_domain()[source]

Check that the mesh covers the FEM computational domain Ω.

Simple necessary (not sufficient) check: whether the convex hull of the mesh vertices covers domain_hull, the convex hull of the union of domain. Does not check the mesh’s actual (possibly non-convex) boundary shape, only its convex extent.

property nbElements

Number of boundary elements (edges, for a 2D mesh).

property nbEdges

Total number of mesh edges.

property geoshape

MFEM geometry type of the mesh elements (e.g. triangle).

property getBoundaryEdge

(nbElements, 2) ndarray of int – vertex-index pairs for each boundary edge.

property boundary_vertex

(nbElements, 2) ndarray of int – vertex-index pairs of the boundary elements (see nbElements).

property fespace

Return the finite elment space

property fespace_order

Vector dimension of fespace (1 for the scalar H1 space used here).

property ndofs

Number of scalar degrees of freedom in fespace (all mesh vertices, before cutting to effective_dofs).

property effective_dofs

Return the effective number of DOFs (associated with vertices)

property GetVSize

Total number of vector DOFs in fespace (ndofs * fespace_order).

property stiff

Sparse stiffness matrix \(G\) (the <grad(phi_i), grad(phi_j)> FEM Laplacian term), shape (effective_dofs, effective_dofs).

property mass

Sparse (diagonal, lumped) mass matrix \(C\), shape (effective_dofs, effective_dofs).

property inner_points

Return the point of the grid

property n_inner_points

Return the point of the grid

property outer_points

Return the boundary point of the grid

property n_outer_points

Return the number of the boundary point of the grid

property totpoints

Return the total point (inner + boundary) of the grid (same as vertex)

property n_totpoints

Return the total point number of the grid

property shape

Tuple (nvertex, nelements, nbElements).

compute_angles(vertex, triangle)[source]

The three interior angles (degrees) of a single triangle.

Parameters:
  • vertex ((n, 2) array_like) – Vertex coordinates array to index into.

  • triangle (array_like of int, length 3) – Indices (into vertex) of the triangle’s three corners.

Returns:

ndarray, shape (3,) – Angles at each of the three corners, in the same order as triangle.

compute_area(vertex, triangle)[source]

Area of a single triangle (shoelace formula).

Parameters:
  • vertex ((n, 2) array_like) – Vertex coordinates array to index into.

  • triangle (array_like of int, length 3) – Indices (into vertex) of the triangle’s three corners.

Returns:

float

compute_stats()[source]

Interior angles and areas of every mesh triangle (via compute_angles/compute_area); used by __init__ (when stats=True) to populate angles/areas and warn about poor-quality (small-angle) meshes.

Returns:

  • angles ((nelements, 3) ndarray) – Interior angles (degrees) of each triangle.

  • areas ((nelements,) ndarray) – Area of each triangle.

property angles

(nelements, 3) ndarray of the triangles’ interior angles (degrees), or None if __init__ was called with stats=False.

property areas

(nelements,) ndarray of triangle areas, or None if __init__ was called with stats=False.

generate_summary()[source]

Build the (label, value) rows shown by summary, as two lists of (str, [str]) pairs (left and right columns).

summary()[source]

Return a statsmodels.iolib.summary.Summary table of mesh/FE diagnostics (vertex/element counts, angle & area ranges, DOFs, matrix shapes, …); also backs __str__.

generate_summary_header()[source]

Cheap mesh-size rows (nelements/nbElements, direct O(1) MFEM accessor calls) - unlike generate_summary(), does not touch nvertex/n_inner_points/n_outer_points/angles/areas, which involve a Python loop over every mesh element or a reduction over the (possibly large) angle/area arrays. Used by __repr__ (via format_info_table).

Returns:

gen_top_left, gen_top_right (list of (str, list)) – Two lists of (label, [value]) rows.

__repr__()[source]

Cheap, side-effect-free representation for interactive use (see generate_summary_header) - see summary()/str() for the full mesh-quality report (angles, areas, inner/outer counts).

print_info(msg)[source]

Print msg prefixed with a UTC timestamp.