"""Module for the discretization class."""
from __future__ import annotations
import abc
from typing import Callable, Type
import numpy as np
import scipy.sparse as sps
import pygeon as pg
[docs]
class Discretization(abc.ABC):
"""
Abstract class for PyGeoN discretization methods.
"""
poly_order: int
"""Polynomial degree of the basis functions"""
tensor_order: int
"""Tensor order of the basis functions"""
[docs]
def __init__(self, keyword: str = pg.UNITARY_DATA) -> None:
"""
Initialize the Discretization object.
Args:
keyword (str): The keyword used to identify the discretization method.
Default is pg.UNITARY_DATA.
Returns:
None
"""
self.keyword = keyword
[docs]
def __repr__(self) -> str:
"""
Return a string representation of the Discretization object.
The string includes the type of the discretization and, if available,
the keyword associated with it.
Returns:
str: A string representation of the Discretization object.
"""
s = f"Discretization of type {self.__class__.__name__}"
if hasattr(self, "keyword"):
s += f" with keyword {self.keyword}"
return s
[docs]
@abc.abstractmethod
def ndof(self, sd: pg.Grid) -> int:
"""
Returns the number of degrees of freedom associated to the method.
Args:
sd: Grid, or a subclass.
Returns:
ndof: the number of degrees of freedom.
"""
[docs]
def assemble_mass_matrix(
self, sd: pg.Grid, data: dict | None = None
) -> sps.csc_array:
"""
Assembles the mass matrix by projecting to the corresponding piecewise
polynomial space.
Args:
sd (pg.Grid): Grid object or a subclass.
data (dict | None): Dictionary with physical parameters for scaling.
Returns:
sps.csc_array: The mass matrix.
"""
return self._assemble_inner_product(sd, data, "assemble_mass_matrix")
[docs]
def assemble_lumped_matrix(
self, sd: pg.Grid, data: dict | None = None
) -> sps.csc_array:
"""
Assembles the lumped mass matrix using the corresponding piecewise polynomial
space.
Args:
sd (pg.Grid): Grid object or a subclass.
data (dict | None): Dictionary with physical parameters for scaling.
Returns:
sps.csc_array: The lumped mass matrix.
"""
return self._assemble_inner_product(sd, data, "assemble_lumped_matrix")
def _assemble_inner_product(
self, sd: pg.Grid, data: dict | None, pwp_method: str
) -> sps.csc_array:
"""
General function to assemble an inner product matrix using the piecewise
polynomial space.
Args:
sd (pg.Grid): Grid object or a subclass.
data (dict | None): Dictionary with physical parameters for scaling.
pwp_method (str): Assembly method of the PiecewisePolynomial space.
Returns:
sps.csc_array: The lumped mass matrix.
"""
# If there are no degrees of freedom, we return an empty matrix.
if self.ndof(sd) == 0:
return sps.csc_array((0, 0))
# We project to the piecewise polynomials and use that mass matrix.
Pi = self.proj_to_PwPolynomials(sd)
pwp = pg.get_PwPolynomials(self.poly_order, self.tensor_order)(self.keyword)
M = getattr(pwp, pwp_method)(sd, data)
return (Pi.T @ M @ Pi).tocsc()
[docs]
@abc.abstractmethod
def assemble_diff_matrix(self, sd: pg.Grid) -> sps.csc_array:
"""
Assembles the matrix corresponding to the differential operator.
Args:
sd (pg.Grid): Grid object or a subclass.
Returns:
sps.csc_array: The differential matrix.
"""
[docs]
def assemble_stiff_matrix(
self, sd: pg.Grid, data: dict | None = None
) -> sps.csc_array:
"""
Assembles the stiffness matrix.
This method takes a grid object `sd` and an optional data dictionary `data` as
input. It first calls the `assemble_diff_matrix` method to obtain the
differential matrix `B`.
Then, it creates an instance of the range discretization class using the `dim`
attribute of `sd`. Next, it calls the `assemble_mass_matrix` method of the range
discretization class to obtain the mass matrix `A`.
Finally, it returns the product of the transpose of `B`, `A`, and `B`.
Args:
sd (pg.Grid): Grid object or a subclass.
data (dict | None): Optional data dictionary. Defaults to None.
Returns:
sps.csc_array: The stiffness matrix.
"""
B = self.assemble_diff_matrix(sd)
discr = self.get_range_discr_class(sd.dim)(self.keyword)
A = discr.assemble_mass_matrix(sd, data)
return (B.T @ A @ B).tocsc()
[docs]
@abc.abstractmethod
def interpolate(
self, sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray]
) -> np.ndarray:
"""
Interpolates a function onto the finite element space
Args:
sd (pg.Grid): Grid, or a subclass.
func (Callable): A function that returns the function values at coordinates.
Returns:
np.ndarray: The values of the degrees of freedom
"""
[docs]
def eval_at_cell_centers(self, sd: pg.Grid) -> sps.csc_array:
"""
Assembles the matrix for evaluating the discretization at the cell centers.
Args:
sd (pg.Grid): Grid object or a subclass.
Returns:
sps.csc_array: The evaluation matrix.
"""
# If there are no degrees of freedom, we return an empty matrix.
if self.ndof(sd) == 0:
return sps.csc_array((pg.AMBIENT_DIM**self.tensor_order, 0))
Pi = self.proj_to_PwPolynomials(sd)
Poly_space = pg.get_PwPolynomials(self.poly_order, self.tensor_order)(
self.keyword
)
return Poly_space.eval_at_cell_centers(sd) @ Pi
[docs]
def source_term(
self, sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray]
) -> np.ndarray:
"""
Assembles the source term by interpolating the given function
and multiplying by the mass matrix
Args:
sd (pg.Grid): Grid object or a subclass.
func (Callable): A function that returns the function values at coordinates.
Returns:
np.ndarray: The evaluation matrix.
"""
return self.assemble_mass_matrix(sd) @ self.interpolate(sd, func)
[docs]
@abc.abstractmethod
def assemble_nat_bc(
self, sd: pg.Grid, func: Callable[[np.ndarray], np.ndarray], b_faces: np.ndarray
) -> np.ndarray:
"""
Assembles the natural boundary condition term
(Tr q, p)_Gamma
Args:
sd (pg.Grid): The grid object.
func (Callable): The function representing the natural boundary condition.
b_faces (np.ndarray): The array of boundary faces.
Returns:
np.ndarray: The assembled natural boundary condition term.
"""
[docs]
@abc.abstractmethod
def get_range_discr_class(self, dim: int) -> Type[pg.Discretization]:
"""
Returns the discretization class that contains the range of the differential
Args:
dim (int): The dimension of the range.
Returns:
pg.Discretization: The discretization class containing the range of the
differential
"""
[docs]
@abc.abstractmethod
def proj_to_PwPolynomials(self, sd: pg.Grid) -> sps.csc_array:
"""
Returns the inclusion matrix that projects the finite element space onto
the lowest order piecewise polynomial space without loss of information.
Args:
sd (pg.Grid): The grid.
Returns:
sps.csc_array: The inclusion matrix.
"""
[docs]
def assemble_broken_grad_matrix(self, sd: pg.Grid) -> sps.csc_array:
"""
Assembles the broken (element-wise) gradient matrix for the given grid.
Args:
sd (pg.Grid): The grid or a subclass.
Returns:
sps.csc_array: The assembled broken gradient matrix.
"""
proj = self.proj_to_PwPolynomials(sd)
pwp = pg.get_PwPolynomials(self.poly_order, self.tensor_order)(self.keyword)
grad = pwp.assemble_broken_grad_matrix(sd)
return grad @ proj
[docs]
def error_l2(
self,
sd: pg.Grid,
num_sol: np.ndarray,
ana_sol: Callable[[np.ndarray], np.ndarray],
relative: bool = True,
poly_order: int | None = -1,
data: dict | None = None,
) -> float:
"""
Returns the l2 error computed against an analytical solution given as a
function. The default behavior interpolates the function as a piecewise
polynomial of order at least 1.
This behavior can be overruled by setting poly_order = None. Then the error is
computed using interpolation, which may result in unexpected superconvergence.
Args:
sd (pg.Grid): Grid, or a subclass.
num_sol (np.ndarray): Vector of the numerical solution.
ana_sol (Callable): Function that represents the analytical solution.
relative (bool, optional): Compute the relative error or not. Defaults to
True.
poly_order (int, optional): If poly_order is specified as an integer,
the error is evaluated in a piecewise polynomial space of that order. If
it is None, we use interpolation.
Returns:
float: The computed error.
"""
# Default case in which we use piecewise linears, at least.
if poly_order == -1:
poly_order = max(self.poly_order, 1)
# If poly_order is None, then we use the interpolant of the space.
if poly_order is None:
int_sol = self.interpolate(sd, ana_sol)
mass = self.assemble_mass_matrix(sd, data)
norm = (int_sol @ mass @ int_sol.T) if relative else 1
diff = num_sol - int_sol
return np.sqrt(diff @ mass @ diff.T / norm)
# If poly_order is specified, we compare the solutions in a piecewise polynomial
# space
else:
poly_space = pg.get_PwPolynomials(poly_order, self.tensor_order)(
self.keyword
)
proj_sol = pg.proj_to_PwPolynomials(self, sd, poly_order) @ num_sol
return poly_space.error_l2(
sd, proj_sol, ana_sol, relative, poly_order=None, data=data
)