From 12c55794453e1e47d390996289230720eeabab5c Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Mon, 8 Jan 2024 12:55:20 -0600 Subject: [PATCH 1/3] Fix ONPolynomialSet normalization --- FIAT/expansions.py | 22 +++++++++++----------- test/unit/test_fiat.py | 25 +++++++++++++++---------- 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index c17cffa0..e944d102 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -9,8 +9,7 @@ import numpy import math -from FIAT import reference_element -from FIAT import jacobi +from FIAT import reference_element, jacobi def morton_index2(p, q=0): @@ -52,7 +51,7 @@ def jacobi_factors(x, y, z, dx, dy, dz): return fa, fb, fc, dfa, dfb, dfc -def dubiner_recurrence(dim, n, order, ref_pts, jacobian): +def dubiner_recurrence(dim, n, order, ref_pts, Jinv, scale): """Dubiner recurrence from (Kirby 2010)""" if order > 2: raise ValueError("Higher order derivatives not supported") @@ -65,8 +64,8 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian): sym_outer = lambda x, y: outer(x, y) + outer(y, x) pad_dim = dim + 2 - dX = pad_jacobian(jacobian, pad_dim) - phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), 1.) + dX = pad_jacobian(Jinv, pad_dim) + phi[0] = sum((ref_pts[i] - ref_pts[i] for i in range(dim)), scale) if dphi is not None: dphi[0] = (phi[0] - phi[0]) * dX[0] if ddphi is not None: @@ -115,9 +114,10 @@ def dubiner_recurrence(dim, n, order, ref_pts, jacobian): c * (fc * ddphi[iprev] + sym_outer(dphi[iprev], dfc) + phi[iprev] * ddfc)) # normalize - for alpha in reference_element.lattice_iter(0, n+1, codim+1): - icur = idx(*alpha) - scale = math.sqrt(sum(alpha) + 0.5 * len(alpha)) + d = codim + 1 + for index in reference_element.lattice_iter(0, n+1, d): + icur = idx(*index) + scale = math.sqrt((2.0 * sum(index) + d) / d) for result in results: result[icur] *= scale return results @@ -165,7 +165,7 @@ def __init__(self, ref_el): v2 = self.base_ref_el.get_vertices() self.A, self.b = reference_element.make_affine_mapping(v1, v2) self.mapping = lambda x: numpy.dot(self.A, x) + self.b - self.scale = numpy.sqrt(numpy.linalg.det(self.A)) + self.scale = numpy.sqrt(1.0 / ref_el.volume()) self._dmats_cache = {} def get_num_members(self, n): @@ -184,7 +184,7 @@ def _tabulate(self, n, pts, order=0): """A version of tabulate() that also works for a single point. """ D = self.ref_el.get_spatial_dimension() - return dubiner_recurrence(D, n, order, self._mapping(pts), self.A) + return dubiner_recurrence(D, n, order, self._mapping(pts), self.A, self.scale) def get_dmats(self, degree): """Returns a numpy array with the expansion coefficients dmat[k, j, i] @@ -289,7 +289,7 @@ def _tabulate(self, n, pts, order=0): vals[i,j] = phi_i(pts[j]), derivs[i,j] = D vals[i,j].""" xs = self._mapping(pts).T results = [] - scale = numpy.sqrt(0.5 + numpy.arange(n+1)) + scale = self.scale * numpy.sqrt(2 * numpy.arange(n+1) + 1) for k in range(order+1): v = numpy.zeros((n + 1, len(xs)), xs.dtype) if n >= k: diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 3a393c6b..323751bf 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -530,16 +530,21 @@ def test_error_point_high_order(element): eval(element) -@pytest.mark.parametrize('cell', [I, T, S]) -def test_expansion_orthonormality(cell): - from FIAT import expansions, quadrature - U = expansions.ExpansionSet(cell) - degree = 10 - rule = quadrature.make_quadrature(cell, degree + 1) - phi = U.tabulate(degree, rule.pts) - w = rule.get_weights() - scale = 0.5 ** -cell.get_spatial_dimension() - results = scale * np.dot(phi, w[:, None] * phi.T) +@pytest.mark.parametrize('degree', [0, 10]) +@pytest.mark.parametrize('dim', range(1, 4)) +@pytest.mark.parametrize('cell', ["default", "ufc"]) +def test_expansion_orthonormality(cell, dim, degree): + from FIAT.expansions import ExpansionSet + from FIAT.quadrature_schemes import create_quadrature + from FIAT import reference_element + make_cell = {"default": reference_element.default_simplex, + "ufc": reference_element.ufc_simplex}[cell] + ref_el = make_cell(dim) + U = ExpansionSet(ref_el) + Q = create_quadrature(ref_el, 2 * degree) + qpts, qwts = Q.get_points(), Q.get_weights() + phi = U.tabulate(degree, qpts) + results = np.dot(np.multiply(phi, qwts), phi.T) assert np.allclose(results, np.eye(results.shape[0])) From 5cec49b66301801f4b71a85f0fb07833d4579848 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 9 Jan 2024 13:47:54 -0600 Subject: [PATCH 2/3] Rescale integral moment dofs to match lowest order point and integral variants --- FIAT/brezzi_douglas_marini.py | 5 ++--- FIAT/expansions.py | 35 ++++++++++++++++++++--------------- FIAT/nedelec.py | 4 ++-- FIAT/nedelec_second_kind.py | 6 +++--- FIAT/polynomial_set.py | 8 ++++---- FIAT/raviart_thomas.py | 7 +++---- 6 files changed, 34 insertions(+), 31 deletions(-) diff --git a/FIAT/brezzi_douglas_marini.py b/FIAT/brezzi_douglas_marini.py index 2c9aeb89..a3f1fbb6 100644 --- a/FIAT/brezzi_douglas_marini.py +++ b/FIAT/brezzi_douglas_marini.py @@ -30,13 +30,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # Facet nodes are \int_F v\cdot n p ds where p \in P_{q} # degree is q Q_ref = create_quadrature(facet, interpolant_deg + degree) - Pq = polynomial_set.ONPolynomialSet(facet, degree) + Pq = polynomial_set.ONPolynomialSet(facet, degree, scale=1) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd - 1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet + n = ref_el.compute_normal(f) phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) diff --git a/FIAT/expansions.py b/FIAT/expansions.py index e944d102..72280922 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -141,31 +141,36 @@ def xi_tetrahedron(eta): class ExpansionSet(object): - def __new__(cls, ref_el, *args, **kwargs): + def __new__(cls, *args, **kwargs): """Returns an ExpansionSet instance appopriate for the given reference element.""" if cls is not ExpansionSet: return super(ExpansionSet, cls).__new__(cls) + ref_el = args[0] if ref_el.get_shape() == reference_element.POINT: - return PointExpansionSet(ref_el) + return PointExpansionSet(*args, **kwargs) elif ref_el.get_shape() == reference_element.LINE: - return LineExpansionSet(ref_el) + return LineExpansionSet(*args, **kwargs) elif ref_el.get_shape() == reference_element.TRIANGLE: - return TriangleExpansionSet(ref_el) + return TriangleExpansionSet(*args, **kwargs) elif ref_el.get_shape() == reference_element.TETRAHEDRON: - return TetrahedronExpansionSet(ref_el) + return TetrahedronExpansionSet(*args, **kwargs) else: raise ValueError("Invalid reference element type.") - def __init__(self, ref_el): + def __init__(self, ref_el, scale=None): + if isinstance(scale, str) and scale.lower() == "l2 piola": + scale = 1.0 / ref_el.volume() + elif scale is None: + scale = math.sqrt(1.0 / ref_el.volume()) self.ref_el = ref_el + self.scale = scale dim = ref_el.get_spatial_dimension() self.base_ref_el = reference_element.default_simplex(dim) v1 = ref_el.get_vertices() v2 = self.base_ref_el.get_vertices() self.A, self.b = reference_element.make_affine_mapping(v1, v2) self.mapping = lambda x: numpy.dot(self.A, x) + self.b - self.scale = numpy.sqrt(1.0 / ref_el.volume()) self._dmats_cache = {} def get_num_members(self, n): @@ -266,10 +271,10 @@ def tabulate_jet(self, n, pts, order=1): class PointExpansionSet(ExpansionSet): """Evaluates the point basis on a point reference element.""" - def __init__(self, ref_el): + def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 0: raise ValueError("Must have a point") - super(PointExpansionSet, self).__init__(ref_el) + super(PointExpansionSet, self).__init__(ref_el, **kwargs) def tabulate(self, n, pts): """Returns a numpy array A[i,j] = phi_i(pts[j]) = 1.0.""" @@ -279,10 +284,10 @@ def tabulate(self, n, pts): class LineExpansionSet(ExpansionSet): """Evaluates the Legendre basis on a line reference element.""" - def __init__(self, ref_el): + def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 1: raise Exception("Must have a line") - super(LineExpansionSet, self).__init__(ref_el) + super(LineExpansionSet, self).__init__(ref_el, **kwargs) def _tabulate(self, n, pts, order=0): """Returns a tuple of (vals, derivs) such that @@ -306,18 +311,18 @@ def _tabulate(self, n, pts, order=0): class TriangleExpansionSet(ExpansionSet): """Evaluates the orthonormal Dubiner basis on a triangular reference element.""" - def __init__(self, ref_el): + def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 2: raise Exception("Must have a triangle") - super(TriangleExpansionSet, self).__init__(ref_el) + super(TriangleExpansionSet, self).__init__(ref_el, **kwargs) class TetrahedronExpansionSet(ExpansionSet): """Collapsed orthonormal polynomial expansion on a tetrahedron.""" - def __init__(self, ref_el): + def __init__(self, ref_el, **kwargs): if ref_el.get_spatial_dimension() != 3: raise Exception("Must be a tetrahedron") - super(TetrahedronExpansionSet, self).__init__(ref_el) + super(TetrahedronExpansionSet, self).__init__(ref_el, **kwargs) def polynomial_dimension(ref_el, degree): diff --git a/FIAT/nedelec.py b/FIAT/nedelec.py index c0513302..c36e2b15 100644 --- a/FIAT/nedelec.py +++ b/FIAT/nedelec.py @@ -123,7 +123,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): if phi_deg >= 0: facet = ref_el.construct_subelement(dim) Q_ref = create_quadrature(facet, interpolant_deg + phi_deg) - Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,)) + Pqmd = polynomial_set.ONPolynomialSet(facet, phi_deg, (dim,), scale="L2 piola") Phis = Pqmd.tabulate(Q_ref.get_points())[(0,) * dim] Phis = numpy.transpose(Phis, (0, 2, 1)) @@ -165,7 +165,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): interpolant_deg = degree cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + phi_deg) - Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg) + Pqmd = polynomial_set.ONPolynomialSet(ref_el, phi_deg, scale="L2 piola") Phis = Pqmd.tabulate(Q.get_points())[(0,) * dim] nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (dim,)) for d in range(dim) for phi in Phis) diff --git a/FIAT/nedelec_second_kind.py b/FIAT/nedelec_second_kind.py index 6e526972..08e5c137 100644 --- a/FIAT/nedelec_second_kind.py +++ b/FIAT/nedelec_second_kind.py @@ -129,7 +129,7 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant ref_facet = cell.construct_subelement(codim) Q_ref = create_quadrature(ref_facet, interpolant_deg + rt_degree) if codim == 1: - Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,)) + Phi = ONPolynomialSet(ref_facet, rt_degree, (codim,), scale="L2 piola") else: # Construct Raviart-Thomas on the reference facet RT = RaviartThomas(ref_facet, rt_degree, variant) @@ -149,10 +149,10 @@ def _generate_facet_dofs(self, codim, cell, degree, offset, variant, interpolant # Get the quadrature and Jacobian on this facet Q_facet = FacetQuadratureRule(cell, codim, facet, Q_ref) J = Q_facet.jacobian() - detJ = Q_facet.jacobian_determinant() + Jdet = Q_facet.jacobian_determinant() # Map Phis -> phis (reference values to physical values) - piola_map = J / detJ + piola_map = J / Jdet phis = numpy.dot(Phis, piola_map.T) phis = numpy.transpose(phis, (0, 2, 1)) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 3a8d376b..9e191a34 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -120,7 +120,7 @@ class ONPolynomialSet(PolynomialSet): """ - def __init__(self, ref_el, degree, shape=tuple()): + def __init__(self, ref_el, degree, shape=tuple(), scale=None): if shape == tuple(): num_components = 1 @@ -130,7 +130,7 @@ def __init__(self, ref_el, degree, shape=tuple()): num_exp_functions = expansions.polynomial_dimension(ref_el, degree) num_members = num_components * num_exp_functions embedded_degree = degree - expansion_set = expansions.ExpansionSet(ref_el) + expansion_set = expansions.ExpansionSet(ref_el, scale=scale) # set up coefficients if shape == tuple(): @@ -211,7 +211,7 @@ class ONSymTensorPolynomialSet(PolynomialSet): """ - def __init__(self, ref_el, degree, size=None): + def __init__(self, ref_el, degree, size=None, scale=None): sd = ref_el.get_spatial_dimension() if size is None: @@ -222,7 +222,7 @@ def __init__(self, ref_el, degree, size=None): num_components = size * (size + 1) // 2 num_members = num_components * num_exp_functions embedded_degree = degree - expansion_set = expansions.ExpansionSet(ref_el) + expansion_set = expansions.ExpansionSet(ref_el, scale=scale) # set up coefficients for symmetric tensors coeffs_shape = (num_members, *shape, num_exp_functions) diff --git a/FIAT/raviart_thomas.py b/FIAT/raviart_thomas.py index c66c4ff4..1f190ce5 100644 --- a/FIAT/raviart_thomas.py +++ b/FIAT/raviart_thomas.py @@ -75,13 +75,12 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1} # degree is q - 1 Q_ref = create_quadrature(facet, interpolant_deg + degree - 1) - Pq = polynomial_set.ONPolynomialSet(facet, degree - 1) + Pq = polynomial_set.ONPolynomialSet(facet, degree - 1, scale=1) Pq_at_qpts = Pq.tabulate(Q_ref.get_points())[(0,)*(sd - 1)] for f in top[sd - 1]: cur = len(nodes) Q = FacetQuadratureRule(ref_el, sd-1, f, Q_ref) - Jdet = Q.jacobian_determinant() - n = ref_el.compute_scaled_normal(f) / Jdet + n = ref_el.compute_normal(f) phis = n[None, :, None] * Pq_at_qpts[:, None, :] nodes.extend(functional.FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis) @@ -91,7 +90,7 @@ def __init__(self, ref_el, degree, variant, interpolant_deg): if degree > 1: cur = len(nodes) Q = create_quadrature(ref_el, interpolant_deg + degree - 2) - Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2) + Pkm1 = polynomial_set.ONPolynomialSet(ref_el, degree - 2, scale=1) Pkm1_at_qpts = Pkm1.tabulate(Q.get_points())[(0,) * sd] nodes.extend(functional.IntegralMoment(ref_el, Q, phi, (d,), (sd,)) for d in range(sd) From ce87b2dde259579830a8e279acb8513c94077392 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 7 Feb 2024 17:10:46 +0000 Subject: [PATCH 3/3] Update FIAT/polynomial_set.py --- FIAT/polynomial_set.py | 1 - 1 file changed, 1 deletion(-) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index eb3edf0b..dd97c486 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -121,7 +121,6 @@ class ONPolynomialSet(PolynomialSet): """ - def __init__(self, ref_el, degree, shape=tuple(), scale=None, variant=None): if shape == tuple(): num_components = 1