"""
Tools for approximating functions
"""
from typing import List
import mpmath # type: ignore
import numpy as np
from ufl.core.expr import Expr # type: ignore
from mtnlion.formula import Arguments, Formula
from mtnlion.tools.cache import persist_to_npy_file
[docs]class Legendre:
"""
Generate Legendre matrices from Eqs. (3.30) and (3.31) in "Continuum-Scale Lithium-Ion Battery Cell Model in FEniCS"
"""
def __init__(self, num_functions: int):
"""
Initialize with the number of polynomials to approximate with
:param num_functions: Number of polynomials to approximate with
"""
self.num_functions = num_functions
self.M_vec = np.vectorize(self.Mmn) # pylint: disable=invalid-name
self.K_vec = np.vectorize(self.Kmn) # pylint: disable=invalid-name
[docs] @staticmethod
def Mmn(row: int, col: int) -> float: # pylint: disable=invalid-name
"""
Calculate the Mmn matrix entries
:param row: row to calculate
:param col: column to calculate
"""
return float(
mpmath.quad(lambda r: r ** 2 * mpmath.legendre(row, 2 * r - 1) * mpmath.legendre(col, 2 * r - 1), [0, 1])
)
[docs] @staticmethod
def Kmn(row: int, col: int) -> float: # pylint: disable=invalid-name
"""
Calculate the Kmn matrix entries
:param row: row to calculate
:param col: column to calculate
"""
def d_legendre(order: int, radius: float) -> float:
"""
Calculate the derivative of the shifted legendre polynomial.
:param order: Order of the polynomial
:param radius: Point to calculate derivative
"""
return mpmath.diff(lambda x: mpmath.legendre(order, 2 * x - 1), radius)
return float(mpmath.quad(lambda r: r ** 2 * d_legendre(row, r) * d_legendre(col, r), [0, 1]))
@property # type: ignore
@persist_to_npy_file(
"legendre_m.npy", lambda cache, cls, *_: cache.shape[0] < cls.num_functions
) # pylint: disable=invalid-name
def M(self):
"""
Calculate the Mmn matrix.
"""
return np.fromfunction(self.M_vec, (self.num_functions, self.num_functions))
@property # type: ignore
@persist_to_npy_file(
"legendre_k.npy", lambda cache, cls, *_: cache.shape[0] < cls.num_functions
) # pylint: disable=invalid-name
def K(self):
"""
Calculate the Kmn matrix
"""
return np.fromfunction(self.K_vec, (self.num_functions, self.num_functions))
[docs]class LagrangeMultiplier(Formula):
"""
This formula provides Lagrange Multiplier functionality for domain boundaries.
"""
def __init__(self, domains: List[str], trial_name: str, lm_name: str = None):
"""
Create a Lagrange Multiplier object that will be evaluated on a given list of domains for a given trial
function.
If the lagrange multiplier name is not provided, it will have the same name as the trial function with "lm_"
prepended.
:param domains: Domains of evaluation
:param trial_name: Name of the trial function this multiplier applies to
:param lm_name: Name of the lagrange multiplier
"""
super(LagrangeMultiplier, self).__init__(domains=domains)
if lm_name is None:
lm_name = f"lm_{trial_name}"
self.name = f"{lm_name}_{self.name}"
self.Variables = self.typedef("Variables", f"{lm_name}, {trial_name}")