Source code for recursivenodes.metrics

'''Additional metrics for evaluating a set of nodes'''

import numpy as np
from recursivenodes.polynomials import (
        proriolkoornwinderdubinervandermonde as pkdv,
        proriolkoornwinderdubinervandermondegrad as pkdvgrad,
        proriolkoornwinderdubinervandermondehess as pkdvhess)
from recursivenodes.utils import coord_map, npolys
from recursivenodes.quadrature import simplexgausslegendre


[docs]def mass_matrix_condition(d, n, nodes, domain='biunit'): '''Compute the condition number of the mass matrix `M_{ij} = \\int_{\\Delta^d} \\varphi_i \\varphi_j\\ dx`, where `\\{\\varphi_i\\}` is the set of shape functions associated with a set of nodes. Args: d (int): The spatial dimension n (int): The polynomial degree nodes (ndarray): The set of polynomial nodes, with shape `(\\binom{n+d}{d}, d)` (unless ``domain=barycentric``, then with shape `(\\binom{n+d}{d}, d+1))` domain (str, optional): The domain of the nodes. See ":ref:`domains`". Returns: float: the condition number `\\kappa_2(M)`. ''' nodes = coord_map(nodes, domain, 'biunit') V = pkdv(d, n, nodes) Vinv = np.linalg.inv(V) M = Vinv.T.dot(Vinv) k = np.linalg.cond(M) return k
[docs]def weak_laplacian_condition(d, n, nodes, domain='biunit'): '''Compute the condition number of the weak-form Laplacian matrix `K_{ij} = \\int_{\\Delta^d} \\nabla \\varphi_i \\cdot \\nabla \\varphi_j\\ dx`, where `\\{\\varphi_i\\}` is the set of shape functions associated with a set of nodes. Args: d (int): The spatial dimension n (int): The polynomial degree nodes (ndarray): The set of polynomial nodes, with shape `(\\binom{n+d}{d}, d)` (unless ``domain=barycentric``, then with shape `(\\binom{n+d}{d}, d+1)`) domain (str, optional): The domain of the nodes. See ":ref:`domains`". The condition number will be computed for the nodes mapped to the biunit simplex. Returns: float: the condition number `\\kappa_2(K)` (with respect to the pseudoinverse). ''' nodes = coord_map(nodes, domain, 'biunit') V = pkdv(d, n, nodes) Vinv = np.linalg.inv(V) q, w = simplexgausslegendre(d, n) npoints = q.ravel().shape[0] // d w = w.reshape((npoints,)) q = q.reshape((npoints, d)) G = pkdvgrad(d, n, q, C=Vinv) K = np.einsum('kil,k,kjl->ij', G, w, G) s = np.linalg.svd(K, compute_uv=False) k = s[0] / s[-2] return k
[docs]def strong_laplacian_condition(d, n, nodes, domain='biunit'): '''Compute the condition number of the strong-form Laplacian matrix `L_{ij} = \\int_{\\Delta^d} \\varphi_i \\cdot \\Delta \\varphi_j\\ dx`, where `\\{\\varphi_i\\}` is the set of shape functions associated with a set of nodes. Args: d (int): The spatial dimension n (int): The polynomial degree nodes (ndarray): The set of polynomial nodes, with shape `(\\binom{n+d}{d}, d)` (unless ``domain=barycentric``, then with shape `(\\binom{n+d}{d}, d+1))` domain (str, optional): The domain of the nodes. See ":ref:`domains`". The condition number will be computed for the nodes mapped to the biunit simplex. Returns: float: the condition number `\\kappa_2(L)` (with respect to the pseudoinverse). ''' nodes = coord_map(nodes, domain, 'biunit') V = pkdv(d, n, nodes) Vinv = np.linalg.inv(V) q, w = simplexgausslegendre(d, n) npoints = q.ravel().shape[0] // d w = w.reshape((npoints,)) if n >= 2: nharm = npolys(d, n) - npolys(d, n-2) else: return 0. q = q.reshape((npoints, d)) B = pkdv(d, n, q, C=Vinv) H = pkdvhess(d, n, q, C=Vinv) L = np.einsum('ki,k,kjll->ij', B, w, H) s = np.linalg.svd(L, compute_uv=False) k = s[0] / s[-(nharm+1)] return k
[docs]def nodal_laplacian_condition(d, n, nodes, domain='biunit'): '''Compute the condition number of the nodal strong-form Laplacian matrix `L_{ij} = \\Delta \\varphi_j(x_i)`, where `\\{\\varphi_i\\}` is the set of shape functions associated with a set of nodes. Args: d (int): The spatial dimension n (int): The polynomial degree nodes (ndarray): The set of polynomial nodes, with shape `(\\binom{n+d}{d}, d)` (unless ``domain=barycentric``, then with shape `(\\binom{n+d}{d}, d+1))` domain (str, optional): The domain of the nodes. See ":ref:`domains`". The condition number will be computed for the nodes mapped to the biunit simplex. Returns: float: the condition number `\\kappa_2(L)` (with respect to the pseudoinverse). ''' nodes = coord_map(nodes, domain, 'biunit') V = pkdv(d, n, nodes) Vinv = np.linalg.inv(V) if n >= 2: nharm = npolys(d, n) - npolys(d, n-2) else: return 0. H = pkdvhess(d, n, nodes, C=Vinv) L = np.einsum('ijll->ij', H) s = np.linalg.svd(L, compute_uv=False) k = s[0] / s[-(nharm+1)] return k
[docs]def nodal_gradient_condition(d, n, nodes, domain='biunit'): '''Compute the condition number of the nodal gradient matrix `G_{(d i + k)j} = \\nabla_k \\varphi_j(x_i)`, where `\\{\\varphi_i\\}` is the set of shape functions associated with a set of nodes. Args: d (int): The spatial dimension n (int): The polynomial degree nodes (ndarray): The set of polynomial nodes, with shape `(\\binom{n+d}{d}, d)` (unless ``domain=barycentric``, then with shape `(\\binom{n+d}{d}, d+1))` domain (str, optional): The domain of the nodes. See ":ref:`domains`". The condition number will be computed for the nodes mapped to the biunit simplex. Returns: float: the condition number `\\kappa_2(G)` (with respect to the pseudoinverse). ''' nodes = coord_map(nodes, domain, 'biunit') V = pkdv(d, n, nodes) Vinv = np.linalg.inv(V) npoints = nodes.shape[0] G = pkdvgrad(d, n, nodes, C=Vinv) G = np.swapaxes(G, 1, 2).reshape((npoints * d, npoints)) s = np.linalg.svd(G, compute_uv=False) k = s[0] / s[-2] return k