From 2f418ec6316579bc318ba1d4dcb63b058958ec8e Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Mon, 18 Mar 2024 20:41:43 -0400 Subject: [PATCH 01/25] implementation of Hu-Zhang element thus far --- FIAT/__init__.py | 2 + FIAT/functional.py | 28 +++++++++ FIAT/hu_zhang.py | 128 ++++++++++++++++++++++++++++++++++++++ test/unit/test_hz.py | 145 +++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 303 insertions(+) create mode 100644 FIAT/hu_zhang.py create mode 100644 test/unit/test_hz.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index ae68cfd72..5bb82a14a 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -37,6 +37,7 @@ from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson from FIAT.arnold_winther import ArnoldWinther from FIAT.arnold_winther import ArnoldWintherNC +from FIAT.hu_zhang import HuZhang from FIAT.mardal_tai_winther import MardalTaiWinther from FIAT.bubble import Bubble, FacetBubble from FIAT.tensor_product import TensorProductElement @@ -94,6 +95,7 @@ "Hellan-Herrmann-Johnson": HellanHerrmannJohnson, "Conforming Arnold-Winther": ArnoldWinther, "Nonconforming Arnold-Winther": ArnoldWintherNC, + "Hu-Zhang": HuZhang, "Mardal-Tai-Winther": MardalTaiWinther} # List of extra elements diff --git a/FIAT/functional.py b/FIAT/functional.py index 44711ca6d..6c7f2f486 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -433,6 +433,12 @@ def __init__(self, cell, entity, mom_deg, comp_deg): super().__init__(cell, n, t, entity, mom_deg, comp_deg, "IntegralNormalTangentialLegendreMoment") +class IntegralLegendreTangentialTangentialMoment(IntegralLegendreBidirectionalMoment): + """Moment of dot(t, dot(tau, t)) against Legendre on entity.""" + def __init__(self, cell, entity, mom_deg, comp_deg): + t = cell.compute_normalized_edge_tangent(entity) + super().__init__(cell, t, t, entity, mom_deg, comp_deg, + "IntegralTangentialTangentialLegendreMoment") class IntegralMomentOfDivergence(Functional): """Functional representing integral of the divergence of the input @@ -631,6 +637,28 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)] super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment") +#class FrancisIntegralMoment(Functional): +# r""" +# internal nodes are \int_K v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0 +# (cmp. Peter Monk - Finite Element Methods for Maxwell's equations p. 129) +# Note that we don't scale by the area of the facet +# +# :arg ref_el: reference element for which F is a codim-0 entity +# :arg Q: quadrature rule on the face +# :arg P_at_qpts: polynomials evaluated at quad points +# :arg facet: which facet. +# """ +# +# def __init__(self, ref_el, Q, P_at_qpts, facet): +# sd = ref_el.get_spatial_dimension() +# weights = Q.get_weights() +# pt_dict = OrderedDict() +# transform = ref_el.get_entity_transform(sd-1, facet) +# pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) +# for pt, wgt, phi in zip(pts, weights, P_at_qpts): +# pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)] +# super().__init__(ref_el, (sd, ), pt_dict, {}, "FrancisIntegralMoment") + class PointScaledNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py new file mode 100644 index 000000000..e0503ae9f --- /dev/null +++ b/FIAT/hu_zhang.py @@ -0,0 +1,128 @@ +# -*- coding: utf-8 -*- +"""Implementation of the Hu-Zhang finite elements.""" + +# Copyright (C) 2024 by Francis Aznaran (University of Notre Dame) +# +# This file is part of FIAT (https://www.fenicsproject.org) +# +# SPDX-License-Identifier: LGPL-3.0-or-later + + +from FIAT.finite_element import CiarletElement +from FIAT.dual_set import DualSet +from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet +from FIAT.functional import ( + PointwiseInnerProductEvaluation as InnerProduct, + FrobeniusIntegralMoment as FIM, + IntegralMomentOfTensorDivergence, + IntegralLegendreNormalNormalMoment, + IntegralLegendreNormalTangentialMoment, + IntegralLegendreTangentialTangentialMoment, + ) + +from FIAT.quadrature import make_quadrature + +from FIAT.bubble import Bubble + +import numpy + + +class HuZhangDual(DualSet): + def __init__(self, cell, degree): + dofs = [] + dof_ids = {} + dof_ids[0] = {0: [], 1: [], 2: []} + dof_ids[1] = {0: [], 1: [], 2: []} + dof_ids[2] = {0: []} + + dof_cur = 0 + + # vertex dofs + vs = cell.get_vertices() + e1 = numpy.array([1.0, 0.0]) + e2 = numpy.array([0.0, 1.0]) + basis = [(e1, e1), (e1, e2), (e2, e2)] + + dof_cur = 0 + + for entity_id in range(3): + node = tuple(vs[entity_id]) + for (v1, v2) in basis: + dofs.append(InnerProduct(cell, v1, v2, node)) + dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3)) + dof_cur += 3 + + # edge dofs now + # moments of normal . sigma against degree p - 2. + for entity_id in range(3): + #for order in (0, degree - 1): #### NB this should also have been range() back with AW! + for order in range(degree - 1): + dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + degree), + IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + degree)] + # NB, mom_deg should actually be order + degree <= 2*degree, but in AW have 6 = 2*degree + dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(degree - 1))) + dof_cur += 2*(degree - 1) + + # internal dofs + Q = make_quadrature(cell, 2*(degree + 1)) + + e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 + e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 + basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices + + # Copying DOFs of Nedelec of 2nd kind (moments against RT) + qs = Q.get_points() + # Create Lagrange bubble nodal basis + CGbubbles = Bubble(cell, degree) + phi = CGbubbles.get_nodal_basis() + + # Evaluate Lagrange bubble basis at quadrature points + + for (v1, v2) in basis: + v1v2t = numpy.outer(v1, v2) + #phi = [phi[i]*v1v2t for i in len(phi)] + fatqp = numpy.zeros((2, 2, len(Q.pts))) + phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + for k in range(len(Q.pts)): + fatqp[:, :, k] = v1v2t + #temp = phiatqpts[k, :] + #fatqp[:, :, k] = temp.reshape((2, 2)) + phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) + dofs.append(FIM(cell, Q, fatqp)) + dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree + 1))) + dof_cur += 3*(degree + 1) + + #for entity_id in range(3): + # for order in range(1, degree): + # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, degree*2)] + + dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*degree*(degree - 1)/2))) + dof_cur += round(3*degree*(degree - 1)/2) + +# # Constraint dofs +# +# Q = make_quadrature(cell, 5) +# +# onp = ONPolynomialSet(cell, 2, (2,)) +# pts = Q.get_points() +# onpvals = onp.tabulate(pts)[0, 0] +# +# for i in list(range(3, 6)) + list(range(9, 12)): +# dofs.append(IntegralMomentOfTensorDivergence(cell, Q, +# onpvals[i, :, :])) +# +# dof_ids[2][0] += list(range(dof_cur, dof_cur + 6)) + + super(HuZhangDual, self).__init__(dofs, cell, dof_ids) + + +class HuZhang(CiarletElement): + """The definition of the Hu-Zhang element. + """ + def __init__(self, cell, degree): + Ps = ONSymTensorPolynomialSet(cell, degree) + Ls = HuZhangDual(cell, degree) + mapping = "double contravariant piola" + super(HuZhang, self).__init__(Ps, Ls, degree, mapping = mapping) diff --git a/test/unit/test_hz.py b/test/unit/test_hz.py new file mode 100644 index 000000000..a61dc91fd --- /dev/null +++ b/test/unit/test_hz.py @@ -0,0 +1,145 @@ +import numpy as np +from FIAT import ufc_simplex, HuZhang, make_quadrature, expansions + + +def test_dofs(): + line = ufc_simplex(1) + T = ufc_simplex(2) + T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)]) + AW = HuZhang(T, 3) + + # check Kronecker property at vertices + + bases = [[[1, 0], [0, 0]], [[0, 1], [1, 0]], [[0, 0], [0, 1]]] + + vert_vals = AW.tabulate(0, T.vertices)[(0, 0)] + for i in range(3): + for j in range(3): + assert np.allclose(vert_vals[3*i+j, :, :, i], bases[j]) + for k in (1, 2): + assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], np.zeros((2, 2))) + + # check edge moments + Qline = make_quadrature(line, 6) + + linebfs = expansions.LineExpansionSet(line) + linevals = linebfs.tabulate(1, Qline.pts) + + # n, n moments + for ed in range(3): + n = T.compute_scaled_normal(ed) + wts = np.asarray(Qline.wts) + nqpline = len(wts) + + vals = AW.tabulate(0, Qline.pts, (1, ed))[(0, 0)] + nnvals = np.zeros((30, nqpline)) + for i in range(30): + for j in range(len(wts)): + nnvals[i, j] = n @ vals[i, :, :, j] @ n + + nnmoments = np.zeros((30, 2)) + + for bf in range(30): + for k in range(nqpline): + for m in (0, 1): + nnmoments[bf, m] += wts[k] * nnvals[bf, k] * linevals[m, k] + + for bf in range(30): + if bf != AW.dual.entity_ids[1][ed][0] and bf != AW.dual.entity_ids[1][ed][2]: + assert np.allclose(nnmoments[bf, :], np.zeros(2)) + + # n, t moments + for ed in range(3): + n = T.compute_scaled_normal(ed) + t = T.compute_edge_tangent(ed) + wts = np.asarray(Qline.wts) + nqpline = len(wts) + + vals = AW.tabulate(0, Qline.pts, (1, ed))[(0, 0)] + ntvals = np.zeros((30, nqpline)) + for i in range(30): + for j in range(len(wts)): + ntvals[i, j] = n @ vals[i, :, :, j] @ t + + ntmoments = np.zeros((30, 2)) + + for bf in range(30): + for k in range(nqpline): + for m in (0, 1): + ntmoments[bf, m] += wts[k] * ntvals[bf, k] * linevals[m, k] + + for bf in range(30): + if bf != AW.dual.entity_ids[1][ed][1] and bf != AW.dual.entity_ids[1][ed][3]: + assert np.allclose(ntmoments[bf, :], np.zeros(2)) + + # check internal dofs + Q = make_quadrature(T, 6) + qpvals = AW.tabulate(0, Q.pts)[(0, 0)] + const_moms = qpvals @ Q.wts + assert np.allclose(const_moms[:21], np.zeros((21, 2, 2))) + assert np.allclose(const_moms[24:], np.zeros((6, 2, 2))) + assert np.allclose(const_moms[21:24, 0, 0], np.asarray([1, 0, 0])) + assert np.allclose(const_moms[21:24, 0, 1], np.asarray([0, 1, 0])) + assert np.allclose(const_moms[21:24, 1, 0], np.asarray([0, 1, 0])) + assert np.allclose(const_moms[21:24, 1, 1], np.asarray([0, 0, 1])) + + +def frob(a, b): + return a.ravel() @ b.ravel() + + +def test_projection(): + T = ufc_simplex(2) + T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)]) + + AW = HuZhang(T, 3) + + Q = make_quadrature(T, 4) + qpts = np.asarray(Q.pts) + qwts = np.asarray(Q.wts) + nqp = len(Q.wts) + + nbf = 24 + m = np.zeros((nbf, nbf)) + b = np.zeros((24,)) + rhs_vals = np.zeros((2, 2, nqp)) + + bfvals = AW.tabulate(0, qpts)[(0, 0)][:nbf, :, :, :] + + for i in range(nbf): + for j in range(nbf): + for k in range(nqp): + m[i, j] += qwts[k] * frob(bfvals[i, :, :, k], + bfvals[j, :, :, k]) + + assert np.linalg.cond(m) < 1.e12 + + comps = [(0, 0), (0, 1), (0, 0)] + + # loop over monomials up to degree 2 + for deg in range(3): + for jj in range(deg+1): + ii = deg-jj + for comp in comps: + b[:] = 0.0 + # set RHS (symmetrically) to be the monomial in + # the proper component. + rhs_vals[comp] = qpts[:, 0]**ii * qpts[:, 1]**jj + rhs_vals[tuple(reversed(comp))] = rhs_vals[comp] + for i in range(nbf): + for k in range(nqp): + b[i] += qwts[k] * frob(bfvals[i, :, :, k], + rhs_vals[:, :, k]) + x = np.linalg.solve(m, b) + + sol_at_qpts = np.zeros(rhs_vals.shape) + for i in range(nbf): + for k in range(nqp): + sol_at_qpts[:, :, k] += x[i] * bfvals[i, :, :, k] + + diff = sol_at_qpts - rhs_vals + err = 0.0 + for k in range(nqp): + err += qwts[k] * frob(diff[:, :, k], diff[:, :, k]) + + assert np.sqrt(err) < 1.e-12 From d184ff47dfb95099ff232199b1e80506f57a1015 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Tue, 19 Mar 2024 14:59:12 -0400 Subject: [PATCH 02/25] Switch to lowest degree = 3 --- FIAT/hu_zhang.py | 48 +++++++++++++++++++++++++++----------------- test/unit/test_hz.py | 18 ++++++++--------- 2 files changed, 39 insertions(+), 27 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index e0503ae9f..ec3ae5ef5 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -56,15 +56,21 @@ def __init__(self, cell, degree): # moments of normal . sigma against degree p - 2. for entity_id in range(3): #for order in (0, degree - 1): #### NB this should also have been range() back with AW! - for order in range(degree - 1): - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + degree), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + degree)] + #for order in range(degree - 1): + for order in range(2): + #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + degree), + # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + degree)] + dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), + IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] # NB, mom_deg should actually be order + degree <= 2*degree, but in AW have 6 = 2*degree - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(degree - 1))) - dof_cur += 2*(degree - 1) + #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(degree - 1))) + dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 4)) + #dof_cur += 2*(degree - 1) + dof_cur += 4 # internal dofs - Q = make_quadrature(cell, 2*(degree + 1)) + #Q = make_quadrature(cell, 2*(degree + 1)) + Q = make_quadrature(cell, 3) e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 @@ -73,33 +79,39 @@ def __init__(self, cell, degree): # Copying DOFs of Nedelec of 2nd kind (moments against RT) qs = Q.get_points() # Create Lagrange bubble nodal basis - CGbubbles = Bubble(cell, degree) + #CGbubbles = Bubble(cell, degree) + CGbubbles = Bubble(cell, 3) phi = CGbubbles.get_nodal_basis() # Evaluate Lagrange bubble basis at quadrature points + + # Copying AWc rather than AWnc internal DOFs, since latter has 4 nested for loops for (v1, v2) in basis: v1v2t = numpy.outer(v1, v2) - #phi = [phi[i]*v1v2t for i in len(phi)] + #phi_times_matrix = [phi[i]*v1v2t for i in len(phi)] fatqp = numpy.zeros((2, 2, len(Q.pts))) + #phiatqpts = numpy.outer(phi_times_matrix.tabulate(qs)[(0,) * 2], v1v2t) phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) for k in range(len(Q.pts)): - fatqp[:, :, k] = v1v2t - #temp = phiatqpts[k, :] - #fatqp[:, :, k] = temp.reshape((2, 2)) - phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) + #fatqp[:, :, k] = v1v2t + temp = phiatqpts[k, :] + fatqp[:, :, k] = temp.reshape((2, 2)) + #phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) dofs.append(FIM(cell, Q, fatqp)) - dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree + 1))) - dof_cur += 3*(degree + 1) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree + 1))) + dof_ids[2][0] = list(range(dof_cur, dof_cur + 3)) + #dof_cur += 3*(degree + 1) + dof_cur += 3 #for entity_id in range(3): # for order in range(1, degree): # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, degree*2)] - dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*degree*(degree - 1)/2))) - dof_cur += round(3*degree*(degree - 1)/2) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*degree*(degree - 1)/2))) + #dof_cur += round(3*degree*(degree - 1)/2) # # Constraint dofs # diff --git a/test/unit/test_hz.py b/test/unit/test_hz.py index a61dc91fd..951df31d2 100644 --- a/test/unit/test_hz.py +++ b/test/unit/test_hz.py @@ -6,13 +6,13 @@ def test_dofs(): line = ufc_simplex(1) T = ufc_simplex(2) T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)]) - AW = HuZhang(T, 3) + HZ = HuZhang(T, 3) # check Kronecker property at vertices bases = [[[1, 0], [0, 0]], [[0, 1], [1, 0]], [[0, 0], [0, 1]]] - vert_vals = AW.tabulate(0, T.vertices)[(0, 0)] + vert_vals = HZ.tabulate(0, T.vertices)[(0, 0)] for i in range(3): for j in range(3): assert np.allclose(vert_vals[3*i+j, :, :, i], bases[j]) @@ -31,7 +31,7 @@ def test_dofs(): wts = np.asarray(Qline.wts) nqpline = len(wts) - vals = AW.tabulate(0, Qline.pts, (1, ed))[(0, 0)] + vals = HZ.tabulate(0, Qline.pts, (1, ed))[(0, 0)] nnvals = np.zeros((30, nqpline)) for i in range(30): for j in range(len(wts)): @@ -45,7 +45,7 @@ def test_dofs(): nnmoments[bf, m] += wts[k] * nnvals[bf, k] * linevals[m, k] for bf in range(30): - if bf != AW.dual.entity_ids[1][ed][0] and bf != AW.dual.entity_ids[1][ed][2]: + if bf != HZ.dual.entity_ids[1][ed][0] and bf != HZ.dual.entity_ids[1][ed][2]: assert np.allclose(nnmoments[bf, :], np.zeros(2)) # n, t moments @@ -55,7 +55,7 @@ def test_dofs(): wts = np.asarray(Qline.wts) nqpline = len(wts) - vals = AW.tabulate(0, Qline.pts, (1, ed))[(0, 0)] + vals = HZ.tabulate(0, Qline.pts, (1, ed))[(0, 0)] ntvals = np.zeros((30, nqpline)) for i in range(30): for j in range(len(wts)): @@ -69,12 +69,12 @@ def test_dofs(): ntmoments[bf, m] += wts[k] * ntvals[bf, k] * linevals[m, k] for bf in range(30): - if bf != AW.dual.entity_ids[1][ed][1] and bf != AW.dual.entity_ids[1][ed][3]: + if bf != HZ.dual.entity_ids[1][ed][1] and bf != HZ.dual.entity_ids[1][ed][3]: assert np.allclose(ntmoments[bf, :], np.zeros(2)) # check internal dofs Q = make_quadrature(T, 6) - qpvals = AW.tabulate(0, Q.pts)[(0, 0)] + qpvals = HZ.tabulate(0, Q.pts)[(0, 0)] const_moms = qpvals @ Q.wts assert np.allclose(const_moms[:21], np.zeros((21, 2, 2))) assert np.allclose(const_moms[24:], np.zeros((6, 2, 2))) @@ -92,7 +92,7 @@ def test_projection(): T = ufc_simplex(2) T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)]) - AW = HuZhang(T, 3) + HZ = HuZhang(T, 3) Q = make_quadrature(T, 4) qpts = np.asarray(Q.pts) @@ -104,7 +104,7 @@ def test_projection(): b = np.zeros((24,)) rhs_vals = np.zeros((2, 2, nqp)) - bfvals = AW.tabulate(0, qpts)[(0, 0)][:nbf, :, :, :] + bfvals = HZ.tabulate(0, qpts)[(0, 0)][:nbf, :, :, :] for i in range(nbf): for j in range(nbf): From eefd2625c630b9f5ef2aedf8d06e773bfd11f89c Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Wed, 20 Mar 2024 11:46:40 -0400 Subject: [PATCH 03/25] More on lowest degree --- FIAT/hu_zhang.py | 39 +++++++++++++++++++++++++++++---------- 1 file changed, 29 insertions(+), 10 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index ec3ae5ef5..87f2132f1 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -22,7 +22,7 @@ from FIAT.quadrature import make_quadrature -from FIAT.bubble import Bubble +from FIAT.bubble import Bubble, FacetBubble # each of these is for the interior DOFs import numpy @@ -35,7 +35,7 @@ def __init__(self, cell, degree): dof_ids[1] = {0: [], 1: [], 2: []} dof_ids[2] = {0: []} - dof_cur = 0 + #dof_cur = 0 # vertex dofs vs = cell.get_vertices() @@ -70,7 +70,7 @@ def __init__(self, cell, degree): # internal dofs #Q = make_quadrature(cell, 2*(degree + 1)) - Q = make_quadrature(cell, 3) + Q = make_quadrature(cell, 1) # In lowest order case I think integration of the product of 2 cubic tensors e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 @@ -84,7 +84,6 @@ def __init__(self, cell, degree): phi = CGbubbles.get_nodal_basis() # Evaluate Lagrange bubble basis at quadrature points - # Copying AWc rather than AWnc internal DOFs, since latter has 4 nested for loops for (v1, v2) in basis: @@ -93,6 +92,7 @@ def __init__(self, cell, degree): fatqp = numpy.zeros((2, 2, len(Q.pts))) #phiatqpts = numpy.outer(phi_times_matrix.tabulate(qs)[(0,) * 2], v1v2t) phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + #print(len(Q.pts)) for k in range(len(Q.pts)): #fatqp[:, :, k] = v1v2t temp = phiatqpts[k, :] @@ -101,17 +101,33 @@ def __init__(self, cell, degree): #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) dofs.append(FIM(cell, Q, fatqp)) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree + 1))) - dof_ids[2][0] = list(range(dof_cur, dof_cur + 3)) - #dof_cur += 3*(degree + 1) - dof_cur += 3 + #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(degree - 1)*(degree - 2)/2)))) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) + #dof_cur += round(3*(degree - 1)*(degree - 2)/2) + #dof_cur += 3 - #for entity_id in range(3): + for entity_id in range(3): # for order in range(1, degree): - # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, degree*2)] + for order in range(1, 3): + # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*degree)] + dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] + + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree - 1)) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) + #dof_cur += 3*(degree - 1) + #dof_cur += 6 + + # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). + # Could instead do via moments against edge bubbles. + #CGEdgeBubbles = FaceBubble() + #for entity_id in range(3): + + # This counting below can be done here, or above for one type of internal DOF at a time #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*degree*(degree - 1)/2))) + dof_ids[2][0] = list(range(dof_cur, dof_cur + 9)) #dof_cur += round(3*degree*(degree - 1)/2) + dof_cur += 9 # # Constraint dofs # @@ -127,6 +143,9 @@ def __init__(self, cell, degree): # # dof_ids[2][0] += list(range(dof_cur, dof_cur + 6)) + #print(dof_cur) + #print(dof_ids) + super(HuZhangDual, self).__init__(dofs, cell, dof_ids) From f545c529a6150c4147119e44ca3f2cc4d34b66a2 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Tue, 26 Mar 2024 20:08:29 -0400 Subject: [PATCH 04/25] Let degree = p. Makes much easier to read --- FIAT/hu_zhang.py | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 87f2132f1..07a4a1bdb 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -29,6 +29,7 @@ class HuZhangDual(DualSet): def __init__(self, cell, degree): + p = degree # This just makes some code below easier to read dofs = [] dof_ids = {} dof_ids[0] = {0: [], 1: [], 2: []} @@ -53,24 +54,25 @@ def __init__(self, cell, degree): dof_cur += 3 # edge dofs now - # moments of normal . sigma against degree p - 2. + # moments of normal component of sigma against degree p - 2. for entity_id in range(3): - #for order in (0, degree - 1): #### NB this should also have been range() back with AW! - #for order in range(degree - 1): + #for order in (0, p - 1): #### NB this should also have been range() back with AW! + #for order in range(p - 1): for order in range(2): - #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + degree), - # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + degree)] + #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + p), + # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + p)] dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] - # NB, mom_deg should actually be order + degree <= 2*degree, but in AW have 6 = 2*degree - #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(degree - 1))) + # NB, mom_deg should actually be order + p <= 2p, but in AW have 6 = 2p + #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(p - 1))) dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 4)) - #dof_cur += 2*(degree - 1) + #dof_cur += 2*(p - 1) dof_cur += 4 # internal dofs - #Q = make_quadrature(cell, 2*(degree + 1)) - Q = make_quadrature(cell, 1) # In lowest order case I think integration of the product of 2 cubic tensors + #Q = make_quadrature(cell, 2*(p + 1)) + #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is + Q = make_quadrature(cell, 3) # In lowest order case I think integration of the product of 2 cubic tensors e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 @@ -79,7 +81,7 @@ def __init__(self, cell, degree): # Copying DOFs of Nedelec of 2nd kind (moments against RT) qs = Q.get_points() # Create Lagrange bubble nodal basis - #CGbubbles = Bubble(cell, degree) + #CGbubbles = Bubble(cell, p) CGbubbles = Bubble(cell, 3) phi = CGbubbles.get_nodal_basis() @@ -101,20 +103,20 @@ def __init__(self, cell, degree): #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) dofs.append(FIM(cell, Q, fatqp)) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(degree - 1)*(degree - 2)/2)))) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(p - 1)*(p - 2)/2)))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) - #dof_cur += round(3*(degree - 1)*(degree - 2)/2) + #dof_cur += round(3*(p - 1)*(p - 2)/2) #dof_cur += 3 for entity_id in range(3): - # for order in range(1, degree): + # for order in range(1, p): for order in range(1, 3): - # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*degree)] + # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(degree - 1)) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1)) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) - #dof_cur += 3*(degree - 1) + #dof_cur += 3*(p - 1) #dof_cur += 6 # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). @@ -124,9 +126,9 @@ def __init__(self, cell, degree): # This counting below can be done here, or above for one type of internal DOF at a time - #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*degree*(degree - 1)/2))) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*p*(p - 1)/2))) dof_ids[2][0] = list(range(dof_cur, dof_cur + 9)) - #dof_cur += round(3*degree*(degree - 1)/2) + #dof_cur += round(3*p*(p - 1)/2) dof_cur += 9 # # Constraint dofs From 89181e4663a3bc145b9146e2225ac6e008bab5a2 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Sun, 26 May 2024 17:33:51 +0100 Subject: [PATCH 05/25] =?UTF-8?q?Nothing=20but=20replace=203=20with=20p?= =?UTF-8?q?=C2=A3?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- FIAT/hu_zhang.py | 48 +++++++++++++++++++++++------------------------- 1 file changed, 23 insertions(+), 25 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 07a4a1bdb..1cc394677 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -26,7 +26,6 @@ import numpy - class HuZhangDual(DualSet): def __init__(self, cell, degree): p = degree # This just makes some code below easier to read @@ -56,23 +55,23 @@ def __init__(self, cell, degree): # edge dofs now # moments of normal component of sigma against degree p - 2. for entity_id in range(3): - #for order in (0, p - 1): #### NB this should also have been range() back with AW! - #for order in range(p - 1): - for order in range(2): - #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + p), - # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + p)] - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] + #for order in (0, p - 1): + for order in range(p - 1): #### NB this should also have been range() back with AW! + #for order in range(2): + dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + p), + IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + p)] + #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), + # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] # NB, mom_deg should actually be order + p <= 2p, but in AW have 6 = 2p - #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(p - 1))) - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 4)) - #dof_cur += 2*(p - 1) - dof_cur += 4 + dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(p - 1))) + #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 4)) + dof_cur += 2*(p - 1) + #dof_cur += 4 # internal dofs #Q = make_quadrature(cell, 2*(p + 1)) - #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is - Q = make_quadrature(cell, 3) # In lowest order case I think integration of the product of 2 cubic tensors + Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is + #Q = make_quadrature(cell, 3) # In lowest order case I think integration of the product of 2 cubic tensors e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 @@ -81,8 +80,8 @@ def __init__(self, cell, degree): # Copying DOFs of Nedelec of 2nd kind (moments against RT) qs = Q.get_points() # Create Lagrange bubble nodal basis - #CGbubbles = Bubble(cell, p) - CGbubbles = Bubble(cell, 3) + CGbubbles = Bubble(cell, p) + #CGbubbles = Bubble(cell, 3) phi = CGbubbles.get_nodal_basis() # Evaluate Lagrange bubble basis at quadrature points @@ -109,10 +108,10 @@ def __init__(self, cell, degree): #dof_cur += 3 for entity_id in range(3): - # for order in range(1, p): - for order in range(1, 3): - # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] - dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] + for order in range(1, p): + # for order in range(1, 3): + dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] + # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1)) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) @@ -124,12 +123,11 @@ def __init__(self, cell, degree): #CGEdgeBubbles = FaceBubble() #for entity_id in range(3): - # This counting below can be done here, or above for one type of internal DOF at a time - #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*p*(p - 1)/2))) - dof_ids[2][0] = list(range(dof_cur, dof_cur + 9)) - #dof_cur += round(3*p*(p - 1)/2) - dof_cur += 9 + dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*p*(p - 1)/2))) + #dof_ids[2][0] = list(range(dof_cur, dof_cur + 9)) + dof_cur += round(3*p*(p - 1)/2) + #dof_cur += 9 # # Constraint dofs # From b0f45c6cd5cf2fed0ee571a04fe4f88f015069f9 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Mon, 27 May 2024 13:20:01 +0100 Subject: [PATCH 06/25] A higher-degree space produced --- FIAT/hu_zhang.py | 49 ++++++++++++++++++++++++++++++++++++------------ 1 file changed, 37 insertions(+), 12 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 1cc394677..02124aec9 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -89,37 +89,62 @@ def __init__(self, cell, degree): for (v1, v2) in basis: v1v2t = numpy.outer(v1, v2) - #phi_times_matrix = [phi[i]*v1v2t for i in len(phi)] - fatqp = numpy.zeros((2, 2, len(Q.pts))) + #phi_times_matrix = [phi[i]*v1v2t for i in range(phi.get_num_members())] + #fatqp = numpy.zeros((2, 2, len(Q.pts))) + Fatqp = numpy.zeros((2, 2, phi.get_num_members())) #phiatqpts = numpy.outer(phi_times_matrix.tabulate(qs)[(0,) * 2], v1v2t) phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) #print(len(Q.pts)) - for k in range(len(Q.pts)): + dim_of_bubbles = phi.get_num_members() + for j in range(dim_of_bubbles): + fatQP = numpy.zeros((2, 2, len(Q.pts))) + # Each DOF here is somehow independent of len(Q.pts) + num_q_pts = len(Q.pts) + for k in range(num_q_pts): + #fatQP[:, :, k] = phiatqpts[j*k:(j + 1)*k, :] + #temp = phiatqpts[j*dim_of_bubbles:(j + 1)*dim_of_bubbles, :] + #temp = phiatqpts[j*k, :] + temp = phiatqpts[j*num_q_pts + k, :] + #print("note: ", temp.shape) + fatQP[:, :, k] = temp.reshape((2, 2)) #fatqp[:, :, k] = v1v2t - temp = phiatqpts[k, :] - fatqp[:, :, k] = temp.reshape((2, 2)) + # temp = phiatqpts[k, :] + #fatqp[:, :, k] = temp.reshape((2, 2)) #phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + dofs.append(FIM(cell, Q, fatQP)) + phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) + #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 0], v1v2t) + #print((phi.tabulate(qs)[(0,) * 2]).shape) + #print(phi_at_qs.shape) + #print(len(phi_at_qs)) + #print(len(phi_at_qs[0, :])) + #print(len(Q.pts)) + #print(phi.get_num_members()) #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) - dofs.append(FIM(cell, Q, fatqp)) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(p - 1)*(p - 2)/2)))) + #dofs.append([FIM(cell, Q, ) for i in range(phi.get_num_members())]) + #Temp = temp.reshape((2, 2)) + #dofs.append(FIM(cell, Q, Temp)) + #dofs.append(FIM(cell, qs, Temp)) + #print(len(FIM(cell, Q, fatqp))) + #dofs.append(FIM(cell, Q, fatqp)) + dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(p - 1)*(p - 2)/2))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += round(3*(p - 1)*(p - 2)/2) #dof_cur += 3 + # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). for entity_id in range(3): for order in range(1, p): # for order in range(1, 3): dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1)) + dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += 3*(p - 1) #dof_cur += 6 - # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). - # Could instead do via moments against edge bubbles. + # Could instead do the tt^T internal dofs via moments against edge bubbles. #CGEdgeBubbles = FaceBubble() #for entity_id in range(3): @@ -145,10 +170,10 @@ def __init__(self, cell, degree): #print(dof_cur) #print(dof_ids) + #print(len(dofs)) super(HuZhangDual, self).__init__(dofs, cell, dof_ids) - class HuZhang(CiarletElement): """The definition of the Hu-Zhang element. """ From 279b9b9a8dd6b855910921f86a27f3ed92849296 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Mon, 27 May 2024 19:07:22 +0100 Subject: [PATCH 07/25] Subtly different internal dof: moment against tt^T, NOT t-t moment. --- FIAT/hu_zhang.py | 35 ++++++++++++++++++++++++++++++++--- 1 file changed, 32 insertions(+), 3 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 02124aec9..3824d3725 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -11,6 +11,7 @@ from FIAT.finite_element import CiarletElement from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet +from FIAT import polynomial_set from FIAT.functional import ( PointwiseInnerProductEvaluation as InnerProduct, FrobeniusIntegralMoment as FIM, @@ -104,7 +105,9 @@ def __init__(self, cell, degree): #fatQP[:, :, k] = phiatqpts[j*k:(j + 1)*k, :] #temp = phiatqpts[j*dim_of_bubbles:(j + 1)*dim_of_bubbles, :] #temp = phiatqpts[j*k, :] - temp = phiatqpts[j*num_q_pts + k, :] + #temp = phiatqpts[j*num_q_pts + k, :] + # NOTE depends how entries of phiatqpts are ordered + temp = phiatqpts[k*dim_of_bubbles + j, :] #print("note: ", temp.shape) fatQP[:, :, k] = temp.reshape((2, 2)) #fatqp[:, :, k] = v1v2t @@ -133,12 +136,36 @@ def __init__(self, cell, degree): #dof_cur += 3 # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). + #ts = cell.compute_tangents( + # Copying BDM + #facet = cell.get_facet_element() + #Q_ref = create_quadrature(facet, p) + Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is + qs = Q.get_points() + #Pp = polynomial_set.ONPolynomialSet(facet, p) + Pp = polynomial_set.ONPolynomialSet(cell, p) + Pp_at_qpts = Pp.tabulate(qs)[(0,) * 2] + #print(Pp_at_qpts.shape) + dim_of_Pp = Pp.get_num_members() # i.e. don't have to call get_nodal_basis() on Pp and then call get_num_members() on that + #print(dim_of_Pp) for entity_id in range(3): + t = cell.compute_edge_tangent(entity_id) + ttT = numpy.outer(t, t) + test_fns_at_qpts = numpy.outer(Pp_at_qpts, ttT) for order in range(1, p): # for order in range(1, 3): - dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] + ## MISTAKE Frobenius inner product with tt^T is not the same as tangential-tangential moment + #dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] - + fatQP = numpy.zeros((2, 2, len(Q.pts))) + num_q_pts = len(Q.pts) + for k in range(num_q_pts): + #temp = phiatqpts[j*num_q_pts + k, :] + # NOTE depends how entries of phiatqpts are ordered, exactly in analogy to the other internal DOFs + temp = test_fns_at_qpts[k*dim_of_bubbles + order, :] + #print("note: ", temp.shape) + fatQP[:, :, k] = temp.reshape((2, 2)) + dofs.append(FIM(cell, Q, fatQP)) dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += 3*(p - 1) @@ -147,6 +174,8 @@ def __init__(self, cell, degree): # Could instead do the tt^T internal dofs via moments against edge bubbles. #CGEdgeBubbles = FaceBubble() #for entity_id in range(3): + + ### NEW complex-based DOFs: airy of fns from the left, div against fns on the right # This counting below can be done here, or above for one type of internal DOF at a time dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*p*(p - 1)/2))) From 917d1360ec09e40226f0e9a69fbbfa778496eab0 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Tue, 11 Jun 2024 16:40:17 -0400 Subject: [PATCH 08/25] Alternative t-t edge dofs --- FIAT/hu_zhang.py | 48 ++++++++++++++++++++++++++++++------------------ 1 file changed, 30 insertions(+), 18 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 3824d3725..5f1f2e4a6 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -12,6 +12,7 @@ from FIAT.dual_set import DualSet from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet from FIAT import polynomial_set +from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import ( PointwiseInnerProductEvaluation as InnerProduct, FrobeniusIntegralMoment as FIM, @@ -139,33 +140,44 @@ def __init__(self, cell, degree): #ts = cell.compute_tangents( # Copying BDM #facet = cell.get_facet_element() - #Q_ref = create_quadrature(facet, p) - Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is - qs = Q.get_points() + #Q = create_quadrature(facet, p) + #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is + #qs = Q.get_points() #Pp = polynomial_set.ONPolynomialSet(facet, p) - Pp = polynomial_set.ONPolynomialSet(cell, p) - Pp_at_qpts = Pp.tabulate(qs)[(0,) * 2] + #Pp = polynomial_set.ONPolynomialSet(cell, p) + #Pp_at_qpts = Pp.tabulate(qs)[(0,) * 2] #print(Pp_at_qpts.shape) - dim_of_Pp = Pp.get_num_members() # i.e. don't have to call get_nodal_basis() on Pp and then call get_num_members() on that + #dim_of_Pp = Pp.get_num_members() # i.e. don't have to call get_nodal_basis() on Pp and then call get_num_members() on that #print(dim_of_Pp) for entity_id in range(3): + pts = cell.make_points(1, entity_id, p + 2) # Gives p + 1 points + print(len(pts), "hi") t = cell.compute_edge_tangent(entity_id) - ttT = numpy.outer(t, t) - test_fns_at_qpts = numpy.outer(Pp_at_qpts, ttT) - for order in range(1, p): + #dofs += [InnerProduct(cell, t, t, pt) for pt in pts] + #ttT = numpy.outer(t, t) + #Q = create_quadrature(entity_id, 2*p - 1) # Since this should give (p - 1) points + #qs = Q.get_points() + #test_fns_at_qpts = numpy.outer(Pp_at_qpts, ttT) + #num_evaluation_pts = len(qs) + #for order in range(1, p): # for order in range(1, 3): ## MISTAKE Frobenius inner product with tt^T is not the same as tangential-tangential moment #dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] - fatQP = numpy.zeros((2, 2, len(Q.pts))) - num_q_pts = len(Q.pts) - for k in range(num_q_pts): - #temp = phiatqpts[j*num_q_pts + k, :] - # NOTE depends how entries of phiatqpts are ordered, exactly in analogy to the other internal DOFs - temp = test_fns_at_qpts[k*dim_of_bubbles + order, :] - #print("note: ", temp.shape) - fatQP[:, :, k] = temp.reshape((2, 2)) - dofs.append(FIM(cell, Q, fatQP)) + #fatQP = numpy.zeros((2, 2, len(Q.pts))) + #num_q_pts = len(Q.pts) + #for k in range(num_q_pts): + # temp = test_fns_at_qpts[order*num_q_pts + k, :] + # # NOTE depends how entries of phiatqpts are ordered, exactly in analogy to the other internal DOFs + # #temp = test_fns_at_qpts[k*dim_of_bubbles + order, :] + # #print("note: ", temp.shape) + # fatQP[:, :, k] = temp.reshape((2, 2)) + #dofs.append(FIM(cell, Q, fatQP)) + P = 0 + for i in range(1, p): + dofs.append(InnerProduct(cell, t, t, pts[i])) + P += 1 + print(P, "here") dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += 3*(p - 1) From 2bf2072208e44c0e851e96f3034f6c78f6198b5b Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Wed, 12 Jun 2024 10:46:46 -0400 Subject: [PATCH 09/25] Temporary commit while some form of HZ appears to be converging at the correct order wrt h --- FIAT/hu_zhang.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 5f1f2e4a6..58f3d8835 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -177,6 +177,8 @@ def __init__(self, cell, degree): for i in range(1, p): dofs.append(InnerProduct(cell, t, t, pts[i])) P += 1 + #for J in range(p + 1): + # print(pts[i]) print(P, "here") dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) From 753a83d615129a2534acc84b9669febb7241e4e8 Mon Sep 17 00:00:00 2001 From: Francis Aznaran Date: Wed, 12 Jun 2024 14:33:51 -0400 Subject: [PATCH 10/25] Alternative interior-bubble DOFs: evaluation at interior nodes --- FIAT/hu_zhang.py | 32 +++++++++++++++++++++++--------- 1 file changed, 23 insertions(+), 9 deletions(-) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 58f3d8835..a4f4db864 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -72,8 +72,9 @@ def __init__(self, cell, degree): # internal dofs #Q = make_quadrature(cell, 2*(p + 1)) - Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is + #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is #Q = make_quadrature(cell, 3) # In lowest order case I think integration of the product of 2 cubic tensors + Q = create_quadrature(cell, 2*p) e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 @@ -96,7 +97,7 @@ def __init__(self, cell, degree): Fatqp = numpy.zeros((2, 2, phi.get_num_members())) #phiatqpts = numpy.outer(phi_times_matrix.tabulate(qs)[(0,) * 2], v1v2t) phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - #print(len(Q.pts)) + #print("length of Q.pts", len(Q.pts)) dim_of_bubbles = phi.get_num_members() for j in range(dim_of_bubbles): fatQP = numpy.zeros((2, 2, len(Q.pts))) @@ -115,15 +116,14 @@ def __init__(self, cell, degree): # temp = phiatqpts[k, :] #fatqp[:, :, k] = temp.reshape((2, 2)) #phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - dofs.append(FIM(cell, Q, fatQP)) + #dofs.append(FIM(cell, Q, fatQP)) phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 0], v1v2t) #print((phi.tabulate(qs)[(0,) * 2]).shape) #print(phi_at_qs.shape) #print(len(phi_at_qs)) #print(len(phi_at_qs[0, :])) - #print(len(Q.pts)) - #print(phi.get_num_members()) + #print(phi.get_num_members(), "bubbles") #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) #dofs.append([FIM(cell, Q, ) for i in range(phi.get_num_members())]) #Temp = temp.reshape((2, 2)) @@ -131,6 +131,17 @@ def __init__(self, cell, degree): #dofs.append(FIM(cell, qs, Temp)) #print(len(FIM(cell, Q, fatqp))) #dofs.append(FIM(cell, Q, fatqp)) + + # Alternative bubble-internal DOFs: just evaluate at the interior points + interior_points = cell.make_points(2, 0, p) # I presume the only cell has entity_id = 0 + num_interior = len(interior_points) + #print("Num interior =", num_interior) + for K in range(num_interior): + #print(interior_points[K]) + for (v1, v2) in basis: + #v1v2t = numpy.outer(v1, v2) + dofs.append(InnerProduct(cell, v1, v2, interior_points[K])) + dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(p - 1)*(p - 2)/2))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += round(3*(p - 1)*(p - 2)/2) @@ -150,8 +161,9 @@ def __init__(self, cell, degree): #dim_of_Pp = Pp.get_num_members() # i.e. don't have to call get_nodal_basis() on Pp and then call get_num_members() on that #print(dim_of_Pp) for entity_id in range(3): - pts = cell.make_points(1, entity_id, p + 2) # Gives p + 1 points - print(len(pts), "hi") + pts = cell.make_points(1, entity_id, p + 2) # Gives p + 1 points. Strange that have to add 2 to the degree? Since p + 1 points determines P_p, not p_{p + 2}? + #pts = cell.make_points(1, entity_id, p) # Could just take p - 1 points, which arises from passing p here + #print(len(pts), "hi") t = cell.compute_edge_tangent(entity_id) #dofs += [InnerProduct(cell, t, t, pt) for pt in pts] #ttT = numpy.outer(t, t) @@ -175,11 +187,13 @@ def __init__(self, cell, degree): #dofs.append(FIM(cell, Q, fatQP)) P = 0 for i in range(1, p): + #for i in range(len(pts)): dofs.append(InnerProduct(cell, t, t, pts[i])) P += 1 #for J in range(p + 1): - # print(pts[i]) - print(P, "here") + # print(pts[J]) + #print(" ") + #print(P, "here") dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) #dof_cur += 3*(p - 1) From f8e12f9aa30bae68e9a45033b396b8483aa37c53 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 8 Oct 2024 16:42:17 +0100 Subject: [PATCH 11/25] Cleanup Hu-Zhang --- FIAT/arnold_winther.py | 8 +- FIAT/functional.py | 28 +--- FIAT/hu_zhang.py | 296 +++++++++++------------------------------ test/unit/test_fiat.py | 7 + test/unit/test_hz.py | 75 +++++------ 5 files changed, 127 insertions(+), 287 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index 5c12fcc76..78a81412e 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -117,17 +117,15 @@ def __init__(self, cell, degree=3): # moments of normal . sigma against constants and linears. for entity_id in range(3): for order in (0, 1): - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] + dofs.append(IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6)) + dofs.append(IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)) + # NB, mom_deg should actually be k + degree <= 2 degree dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) dof_cur += 4 # internal dofs: constant moments of three unique components Q = make_quadrature(cell, 3) - e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 - e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 - basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices for (v1, v2) in basis: v1v2t = numpy.outer(v1, v2) fatqp = numpy.zeros((2, 2, len(Q.pts))) diff --git a/FIAT/functional.py b/FIAT/functional.py index af1df4964..2ad24dab7 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -169,7 +169,7 @@ class PointEvaluation(Functional): particular point x.""" def __init__(self, ref_el, x): - pt_dict = {x: [(1.0, tuple())]} + pt_dict = {tuple(x): [(1.0, tuple())]} Functional.__init__(self, ref_el, tuple(), pt_dict, {}, "PointEval") def __call__(self, fn): @@ -414,6 +414,7 @@ def __init__(self, cell, entity, mom_deg, comp_deg): super().__init__(cell, n, t, entity, mom_deg, comp_deg, "IntegralNormalTangentialLegendreMoment") + class IntegralLegendreTangentialTangentialMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(t, dot(tau, t)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): @@ -421,6 +422,7 @@ def __init__(self, cell, entity, mom_deg, comp_deg): super().__init__(cell, t, t, entity, mom_deg, comp_deg, "IntegralTangentialTangentialLegendreMoment") + class IntegralMomentOfDivergence(Functional): """Functional representing integral of the divergence of the input against some tabulated function f.""" @@ -605,28 +607,6 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): for pt, wt in zip(pts, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment") -#class FrancisIntegralMoment(Functional): -# r""" -# internal nodes are \int_K v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0 -# (cmp. Peter Monk - Finite Element Methods for Maxwell's equations p. 129) -# Note that we don't scale by the area of the facet -# -# :arg ref_el: reference element for which F is a codim-0 entity -# :arg Q: quadrature rule on the face -# :arg P_at_qpts: polynomials evaluated at quad points -# :arg facet: which facet. -# """ -# -# def __init__(self, ref_el, Q, P_at_qpts, facet): -# sd = ref_el.get_spatial_dimension() -# weights = Q.get_weights() -# pt_dict = OrderedDict() -# transform = ref_el.get_entity_transform(sd-1, facet) -# pts = tuple(map(lambda p: tuple(transform(p)), Q.get_points())) -# for pt, wgt, phi in zip(pts, weights, P_at_qpts): -# pt_dict[pt] = [(wgt*phi[i], (i, )) for i in range(sd)] -# super().__init__(ref_el, (sd, ), pt_dict, {}, "FrancisIntegralMoment") - class PointScaledNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a @@ -683,7 +663,7 @@ def __init__(self, ref_el, v, w, pt): wvT = numpy.outer(w, v) shp = wvT.shape - pt_dict = {pt: [(wvT[idx], idx) for idx in index_iterator(shp)]} + pt_dict = {tuple(pt): [(wvT[idx], idx) for idx in index_iterator(shp)]} super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval") diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index a4f4db864..8eed4c919 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -8,234 +8,98 @@ # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT.finite_element import CiarletElement -from FIAT.dual_set import DualSet -from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet -from FIAT import polynomial_set +from FIAT import finite_element, polynomial_set, dual_set +from FIAT.reference_element import TRIANGLE +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -from FIAT.functional import ( - PointwiseInnerProductEvaluation as InnerProduct, - FrobeniusIntegralMoment as FIM, - IntegralMomentOfTensorDivergence, - IntegralLegendreNormalNormalMoment, - IntegralLegendreNormalTangentialMoment, - IntegralLegendreTangentialTangentialMoment, - ) - -from FIAT.quadrature import make_quadrature - -from FIAT.bubble import Bubble, FacetBubble # each of these is for the interior DOFs +from FIAT.functional import (PointwiseInnerProductEvaluation as InnerProduct, + FrobeniusIntegralMoment as FIM, + IntegralMomentOfTensorDivergence, + IntegralLegendreNormalNormalMoment, + IntegralLegendreNormalTangentialMoment, + IntegralLegendreTangentialTangentialMoment, + ) import numpy -class HuZhangDual(DualSet): - def __init__(self, cell, degree): - p = degree # This just makes some code below easier to read - dofs = [] - dof_ids = {} - dof_ids[0] = {0: [], 1: [], 2: []} - dof_ids[1] = {0: [], 1: [], 2: []} - dof_ids[2] = {0: []} - #dof_cur = 0 +class HuZhangDual(dual_set.DualSet): + def __init__(self, ref_el, degree, variant): + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] # vertex dofs - vs = cell.get_vertices() e1 = numpy.array([1.0, 0.0]) e2 = numpy.array([0.0, 1.0]) basis = [(e1, e1), (e1, e2), (e2, e2)] - - dof_cur = 0 - - for entity_id in range(3): - node = tuple(vs[entity_id]) - for (v1, v2) in basis: - dofs.append(InnerProduct(cell, v1, v2, node)) - dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3)) - dof_cur += 3 - - # edge dofs now - # moments of normal component of sigma against degree p - 2. - for entity_id in range(3): - #for order in (0, p - 1): - for order in range(p - 1): #### NB this should also have been range() back with AW! - #for order in range(2): - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, order + p), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, order + p)] - #dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), - # IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] - # NB, mom_deg should actually be order + p <= 2p, but in AW have 6 = 2p - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 2*(p - 1))) - #dof_ids[1][entity_id] = list(range(dof_cur, dof_cur + 4)) - dof_cur += 2*(p - 1) - #dof_cur += 4 + for entity_id in sorted(top[0]): + cur = len(nodes) + pt, = ref_el.make_points(0, entity_id, degree) + nodes.extend(InnerProduct(ref_el, v1, v2, pt) for (v1, v2) in basis) + entity_ids[0][entity_id].extend(range(cur, len(nodes))) + + # edge dofs: moments of normal component of sigma against degree p - 2. + for entity_id in sorted(top[1]): + cur = len(nodes) + for k in range(degree-1): + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity_id, k, k + degree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity_id, k, k + degree)) + # NB, mom_deg should actually be k + degree <= 2 degree, but in AW have 6 = 2 degree + entity_ids[1][entity_id].extend(range(cur, len(nodes))) # internal dofs - #Q = make_quadrature(cell, 2*(p + 1)) - #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is - #Q = make_quadrature(cell, 3) # In lowest order case I think integration of the product of 2 cubic tensors - Q = create_quadrature(cell, 2*p) - - e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 - e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 - basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices - - # Copying DOFs of Nedelec of 2nd kind (moments against RT) - qs = Q.get_points() - # Create Lagrange bubble nodal basis - CGbubbles = Bubble(cell, p) - #CGbubbles = Bubble(cell, 3) - phi = CGbubbles.get_nodal_basis() - - # Evaluate Lagrange bubble basis at quadrature points - # Copying AWc rather than AWnc internal DOFs, since latter has 4 nested for loops - - for (v1, v2) in basis: - v1v2t = numpy.outer(v1, v2) - #phi_times_matrix = [phi[i]*v1v2t for i in range(phi.get_num_members())] - #fatqp = numpy.zeros((2, 2, len(Q.pts))) - Fatqp = numpy.zeros((2, 2, phi.get_num_members())) - #phiatqpts = numpy.outer(phi_times_matrix.tabulate(qs)[(0,) * 2], v1v2t) - phiatqpts = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - #print("length of Q.pts", len(Q.pts)) - dim_of_bubbles = phi.get_num_members() - for j in range(dim_of_bubbles): - fatQP = numpy.zeros((2, 2, len(Q.pts))) - # Each DOF here is somehow independent of len(Q.pts) - num_q_pts = len(Q.pts) - for k in range(num_q_pts): - #fatQP[:, :, k] = phiatqpts[j*k:(j + 1)*k, :] - #temp = phiatqpts[j*dim_of_bubbles:(j + 1)*dim_of_bubbles, :] - #temp = phiatqpts[j*k, :] - #temp = phiatqpts[j*num_q_pts + k, :] - # NOTE depends how entries of phiatqpts are ordered - temp = phiatqpts[k*dim_of_bubbles + j, :] - #print("note: ", temp.shape) - fatQP[:, :, k] = temp.reshape((2, 2)) - #fatqp[:, :, k] = v1v2t - # temp = phiatqpts[k, :] - #fatqp[:, :, k] = temp.reshape((2, 2)) - #phi_at_qs[:, :, k] = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - #dofs.append(FIM(cell, Q, fatQP)) - phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 2], v1v2t) - #phi_at_qs = numpy.outer(phi.tabulate(qs)[(0,) * 0], v1v2t) - #print((phi.tabulate(qs)[(0,) * 2]).shape) - #print(phi_at_qs.shape) - #print(len(phi_at_qs)) - #print(len(phi_at_qs[0, :])) - #print(phi.get_num_members(), "bubbles") - #dofs.append([FIM(cell, Q, phi_at_qs[i, :]) for i in range(len(phi_at_qs))]) - #dofs.append([FIM(cell, Q, ) for i in range(phi.get_num_members())]) - #Temp = temp.reshape((2, 2)) - #dofs.append(FIM(cell, Q, Temp)) - #dofs.append(FIM(cell, qs, Temp)) - #print(len(FIM(cell, Q, fatqp))) - #dofs.append(FIM(cell, Q, fatqp)) - - # Alternative bubble-internal DOFs: just evaluate at the interior points - interior_points = cell.make_points(2, 0, p) # I presume the only cell has entity_id = 0 - num_interior = len(interior_points) - #print("Num interior =", num_interior) - for K in range(num_interior): - #print(interior_points[K]) + cur = len(nodes) + if variant == "integral": + Q = create_quadrature(ref_el, 2*degree-sd-1) + # Moments against P_{degree-3} for each component + phi = polynomial_set.ONPolynomialSet(ref_el, degree-sd-1) + phi_at_qpts = phi.tabulate(Q.get_points())[(0,) * sd] for (v1, v2) in basis: - #v1v2t = numpy.outer(v1, v2) - dofs.append(InnerProduct(cell, v1, v2, interior_points[K])) - - dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*(p - 1)*(p - 2)/2))) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) - #dof_cur += round(3*(p - 1)*(p - 2)/2) - #dof_cur += 3 - - # More internal dofs: evaluation of interior-of-edge Lagrange functions, inner product with tt^T for each edge. Note these are evaluated on the edge, but not shared between cells (hence internal). - #ts = cell.compute_tangents( - # Copying BDM - #facet = cell.get_facet_element() - #Q = create_quadrature(facet, p) - #Q = make_quadrature(cell, p) # p points -> exactly integrate polys of degree 2p + 1 -> in particular a product of two degree p things, which is what this DOF is - #qs = Q.get_points() - #Pp = polynomial_set.ONPolynomialSet(facet, p) - #Pp = polynomial_set.ONPolynomialSet(cell, p) - #Pp_at_qpts = Pp.tabulate(qs)[(0,) * 2] - #print(Pp_at_qpts.shape) - #dim_of_Pp = Pp.get_num_members() # i.e. don't have to call get_nodal_basis() on Pp and then call get_num_members() on that - #print(dim_of_Pp) - for entity_id in range(3): - pts = cell.make_points(1, entity_id, p + 2) # Gives p + 1 points. Strange that have to add 2 to the degree? Since p + 1 points determines P_p, not p_{p + 2}? - #pts = cell.make_points(1, entity_id, p) # Could just take p - 1 points, which arises from passing p here - #print(len(pts), "hi") - t = cell.compute_edge_tangent(entity_id) - #dofs += [InnerProduct(cell, t, t, pt) for pt in pts] - #ttT = numpy.outer(t, t) - #Q = create_quadrature(entity_id, 2*p - 1) # Since this should give (p - 1) points - #qs = Q.get_points() - #test_fns_at_qpts = numpy.outer(Pp_at_qpts, ttT) - #num_evaluation_pts = len(qs) - #for order in range(1, p): - # for order in range(1, 3): - ## MISTAKE Frobenius inner product with tt^T is not the same as tangential-tangential moment - #dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 2*p)] - # dofs += [IntegralLegendreTangentialTangentialMoment(cell, entity_id, order, 6)] - #fatQP = numpy.zeros((2, 2, len(Q.pts))) - #num_q_pts = len(Q.pts) - #for k in range(num_q_pts): - # temp = test_fns_at_qpts[order*num_q_pts + k, :] - # # NOTE depends how entries of phiatqpts are ordered, exactly in analogy to the other internal DOFs - # #temp = test_fns_at_qpts[k*dim_of_bubbles + order, :] - # #print("note: ", temp.shape) - # fatQP[:, :, k] = temp.reshape((2, 2)) - #dofs.append(FIM(cell, Q, fatQP)) - P = 0 - for i in range(1, p): - #for i in range(len(pts)): - dofs.append(InnerProduct(cell, t, t, pts[i])) - P += 1 - #for J in range(p + 1): - # print(pts[J]) - #print(" ") - #print(P, "here") - dof_ids[2][0] = list(range(dof_cur, dof_cur + 3*(p - 1))) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 6)) - #dof_cur += 3*(p - 1) - #dof_cur += 6 - - # Could instead do the tt^T internal dofs via moments against edge bubbles. - #CGEdgeBubbles = FaceBubble() - #for entity_id in range(3): - - ### NEW complex-based DOFs: airy of fns from the left, div against fns on the right - - # This counting below can be done here, or above for one type of internal DOF at a time - dof_ids[2][0] = list(range(dof_cur, dof_cur + round(3*p*(p - 1)/2))) - #dof_ids[2][0] = list(range(dof_cur, dof_cur + 9)) - dof_cur += round(3*p*(p - 1)/2) - #dof_cur += 9 - -# # Constraint dofs -# -# Q = make_quadrature(cell, 5) -# -# onp = ONPolynomialSet(cell, 2, (2,)) -# pts = Q.get_points() -# onpvals = onp.tabulate(pts)[0, 0] -# -# for i in list(range(3, 6)) + list(range(9, 12)): -# dofs.append(IntegralMomentOfTensorDivergence(cell, Q, -# onpvals[i, :, :])) -# -# dof_ids[2][0] += list(range(dof_cur, dof_cur + 6)) - - #print(dof_cur) - #print(dof_ids) - #print(len(dofs)) - - super(HuZhangDual, self).__init__(dofs, cell, dof_ids) - -class HuZhang(CiarletElement): + Phi_at_qpts = numpy.outer(v1, v2)[None, :, :, None] * phi_at_qpts[:, None, None, :] + nodes.extend(FIM(ref_el, Q, Phi) for Phi in Phi_at_qpts) + + # More internal dofs: tangential-tangential moments against bubbles for each edge + # Note these are evaluated on the edge, but not shared between cells (hence internal). + facet = ref_el.get_facet_element() + Q = create_quadrature(facet, 2*degree-sd) + phi = polynomial_set.ONPolynomialSet(facet, degree-sd) + phi_at_qpts = phi.tabulate(Q.get_points())[(0,) * (sd-1)] + for entity_id in sorted(top[1]): + Q_mapped = FacetQuadratureRule(ref_el, sd-1, entity_id, Q) + t = ref_el.compute_edge_tangent(entity_id) + Phi_at_qpts = numpy.outer(t, t)[None, :, :, None] * phi_at_qpts[:, None, None, :] + nodes.extend(FIM(ref_el, Q_mapped, Phi) for Phi in Phi_at_qpts) + + elif variant == "point": + # Evaluation at interior points for each component + interior_points = ref_el.make_points(sd, 0, degree) + nodes.extend(InnerProduct(ref_el, v1, v2, pt) + for pt in interior_points for (v1, v2) in basis) + + # More internal dofs: tangential-tangential point evaluations. + # Note these are evaluated on the edge, but not shared between cells (hence internal). + for entity_id in sorted(top[1]): + t = ref_el.compute_edge_tangent(entity_id) + pts = ref_el.make_points(sd-1, entity_id, degree) + nodes.extend(InnerProduct(ref_el, t, t, pt) for pt in pts) + else: + raise ValueError(f"Unsupported variant {variant}") + + entity_ids[2][0].extend(range(cur, len(nodes))) + super().__init__(nodes, ref_el, entity_ids) + + +class HuZhang(finite_element.CiarletElement): """The definition of the Hu-Zhang element. """ - def __init__(self, cell, degree): - Ps = ONSymTensorPolynomialSet(cell, degree) - Ls = HuZhangDual(cell, degree) - mapping = "double contravariant piola" - super(HuZhang, self).__init__(Ps, Ls, degree, mapping = mapping) + def __init__(self, ref_el, degree=3, variant="integral"): + if degree < 3: + raise ValueError("Hu-Zhang only defined for degree >= 3") + if ref_el.shape != TRIANGLE: + raise ValueError("Hu-Zhang only defined on triangles") + poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + dual = HuZhangDual(ref_el, degree, variant) + formdegree = ref_el.get_spatial_dimension() - 1 + super().__init__(poly_set, dual, degree, formdegree=formdegree, mapping="double contravariant piola") diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 5e2cf264d..2f007a562 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -42,6 +42,8 @@ from FIAT.tensor_product import TensorProductElement # noqa: F401 from FIAT.tensor_product import FlattenedDimensions # noqa: F401 from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 +from FIAT.arnold_winther import ArnoldWinther, ArnoldWintherNC # noqa: F401 +from FIAT.hu_zhang import HuZhang # noqa: F401 from FIAT.bernardi_raugel import BernardiRaugel # noqa: F401 from FIAT.argyris import Argyris # noqa: F401 from FIAT.hermite import CubicHermite # noqa: F401 @@ -307,6 +309,11 @@ def __init__(self, a, b): "Morley(T)", "BernardiRaugel(T)", "BernardiRaugel(S)", + "ArnoldWinther(T, 3)", + "HuZhang(T, 3)", + "HuZhang(T, 4)", + "HuZhang(T, 3, 'point')", + "HuZhang(T, 4, 'point')", # Macroelements "Lagrange(T, 1, 'iso')", diff --git a/test/unit/test_hz.py b/test/unit/test_hz.py index 951df31d2..0a761f40e 100644 --- a/test/unit/test_hz.py +++ b/test/unit/test_hz.py @@ -1,5 +1,8 @@ +import pytest import numpy as np -from FIAT import ufc_simplex, HuZhang, make_quadrature, expansions + +from FIAT import ufc_simplex, HuZhang, expansions +from FIAT.quadrature_schemes import create_quadrature def test_dofs(): @@ -17,10 +20,10 @@ def test_dofs(): for j in range(3): assert np.allclose(vert_vals[3*i+j, :, :, i], bases[j]) for k in (1, 2): - assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], np.zeros((2, 2))) + assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], 0) # check edge moments - Qline = make_quadrature(line, 6) + Qline = create_quadrature(line, 6) linebfs = expansions.LineExpansionSet(line) linevals = linebfs.tabulate(1, Qline.pts) @@ -46,7 +49,7 @@ def test_dofs(): for bf in range(30): if bf != HZ.dual.entity_ids[1][ed][0] and bf != HZ.dual.entity_ids[1][ed][2]: - assert np.allclose(nnmoments[bf, :], np.zeros(2)) + assert np.allclose(nnmoments[bf, :], 0) # n, t moments for ed in range(3): @@ -70,76 +73,64 @@ def test_dofs(): for bf in range(30): if bf != HZ.dual.entity_ids[1][ed][1] and bf != HZ.dual.entity_ids[1][ed][3]: - assert np.allclose(ntmoments[bf, :], np.zeros(2)) + assert np.allclose(ntmoments[bf, :], 0) # check internal dofs - Q = make_quadrature(T, 6) + Q = create_quadrature(T, 6) qpvals = HZ.tabulate(0, Q.pts)[(0, 0)] const_moms = qpvals @ Q.wts - assert np.allclose(const_moms[:21], np.zeros((21, 2, 2))) - assert np.allclose(const_moms[24:], np.zeros((6, 2, 2))) - assert np.allclose(const_moms[21:24, 0, 0], np.asarray([1, 0, 0])) - assert np.allclose(const_moms[21:24, 0, 1], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[21:24, 1, 0], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[21:24, 1, 1], np.asarray([0, 0, 1])) + assert np.allclose(const_moms[:21], 0) + assert np.allclose(const_moms[24:], 0) + assert np.allclose(const_moms[21:24, 0, 0], [1, 0, 0]) + assert np.allclose(const_moms[21:24, 0, 1], [0, 1, 0]) + assert np.allclose(const_moms[21:24, 1, 0], [0, 1, 0]) + assert np.allclose(const_moms[21:24, 1, 1], [0, 0, 1]) def frob(a, b): return a.ravel() @ b.ravel() -def test_projection(): +@pytest.mark.parametrize("variant", ("integral", "point")) +def test_projection(variant): T = ufc_simplex(2) T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)]) - HZ = HuZhang(T, 3) + p = 3 + HZ = HuZhang(T, p, variant) - Q = make_quadrature(T, 4) + Q = create_quadrature(T, 6) qpts = np.asarray(Q.pts) qwts = np.asarray(Q.wts) nqp = len(Q.wts) - nbf = 24 - m = np.zeros((nbf, nbf)) - b = np.zeros((24,)) - rhs_vals = np.zeros((2, 2, nqp)) + nbf = HZ.space_dimension() - 3 * (p-1) + rhs_vals = np.zeros((1, 2, 2, nqp)) bfvals = HZ.tabulate(0, qpts)[(0, 0)][:nbf, :, :, :] - - for i in range(nbf): - for j in range(nbf): - for k in range(nqp): - m[i, j] += qwts[k] * frob(bfvals[i, :, :, k], - bfvals[j, :, :, k]) + ells = np.multiply(bfvals, qwts) + m = np.tensordot(ells, bfvals, (range(1, ells.ndim),)*2) assert np.linalg.cond(m) < 1.e12 - comps = [(0, 0), (0, 1), (0, 0)] + comps = [(0, 0), (0, 1), (1, 1)] # loop over monomials up to degree 2 for deg in range(3): for jj in range(deg+1): ii = deg-jj for comp in comps: - b[:] = 0.0 # set RHS (symmetrically) to be the monomial in # the proper component. - rhs_vals[comp] = qpts[:, 0]**ii * qpts[:, 1]**jj - rhs_vals[tuple(reversed(comp))] = rhs_vals[comp] - for i in range(nbf): - for k in range(nqp): - b[i] += qwts[k] * frob(bfvals[i, :, :, k], - rhs_vals[:, :, k]) - x = np.linalg.solve(m, b) + rhs_vals[...] = 0 + rhs_vals[0][comp] = qpts[:, 0]**ii * qpts[:, 1]**jj + rhs_vals[0][tuple(reversed(comp))] = rhs_vals[0][comp] - sol_at_qpts = np.zeros(rhs_vals.shape) - for i in range(nbf): - for k in range(nqp): - sol_at_qpts[:, :, k] += x[i] * bfvals[i, :, :, k] + b = np.tensordot(ells, rhs_vals, (range(1, ells.ndim),)*2) + x = np.linalg.solve(m, b) - diff = sol_at_qpts - rhs_vals - err = 0.0 - for k in range(nqp): - err += qwts[k] * frob(diff[:, :, k], diff[:, :, k]) + sol_at_qpts = np.tensordot(x, bfvals, (0, 0)) + diff = (sol_at_qpts - rhs_vals)**2 + err = np.linalg.norm(np.tensordot(diff, qwts, (-1, -1))[0], "fro") assert np.sqrt(err) < 1.e-12 From 4bd33035786410a755df7edfa079f41eeb337380 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 00:21:34 +0100 Subject: [PATCH 12/25] Fix point variant --- FIAT/arnold_winther.py | 9 ++- FIAT/hu_zhang.py | 64 ++++++++----------- test/unit/test_hz.py | 136 ----------------------------------------- 3 files changed, 30 insertions(+), 179 deletions(-) delete mode 100644 test/unit/test_hz.py diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index 78a81412e..84febd88c 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -41,8 +41,8 @@ def __init__(self, cell, degree=2): # moments of normal . sigma against constants and linears. for entity_id in range(3): # a triangle has 3 edges for order in (0, 1): - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] + dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, degree), + IntegralLegendreNormalTangentialMoment(cell, entity_id, order, degree)] dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) dof_cur += 4 @@ -117,9 +117,8 @@ def __init__(self, cell, degree=3): # moments of normal . sigma against constants and linears. for entity_id in range(3): for order in (0, 1): - dofs.append(IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6)) - dofs.append(IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)) - # NB, mom_deg should actually be k + degree <= 2 degree + dofs.append(IntegralLegendreNormalNormalMoment(cell, entity_id, order, degree)) + dofs.append(IntegralLegendreNormalTangentialMoment(cell, entity_id, order, degree)) dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) dof_cur += 4 diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 8eed4c919..83b486116 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -9,22 +9,22 @@ from FIAT import finite_element, polynomial_set, dual_set +from FIAT.check_format_variant import check_format_variant from FIAT.reference_element import TRIANGLE -from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (PointwiseInnerProductEvaluation as InnerProduct, FrobeniusIntegralMoment as FIM, - IntegralMomentOfTensorDivergence, IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment, - IntegralLegendreTangentialTangentialMoment, ) import numpy class HuZhangDual(dual_set.DualSet): - def __init__(self, ref_el, degree, variant): + def __init__(self, ref_el, degree, variant, qdegree): + if qdegree is None: + qdegree = degree top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} @@ -44,46 +44,33 @@ def __init__(self, ref_el, degree, variant): for entity_id in sorted(top[1]): cur = len(nodes) for k in range(degree-1): - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity_id, k, k + degree)) - nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity_id, k, k + degree)) - # NB, mom_deg should actually be k + degree <= 2 degree, but in AW have 6 = 2 degree + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity_id, k, qdegree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity_id, k, qdegree)) entity_ids[1][entity_id].extend(range(cur, len(nodes))) # internal dofs cur = len(nodes) - if variant == "integral": - Q = create_quadrature(ref_el, 2*degree-sd-1) - # Moments against P_{degree-3} for each component - phi = polynomial_set.ONPolynomialSet(ref_el, degree-sd-1) - phi_at_qpts = phi.tabulate(Q.get_points())[(0,) * sd] - for (v1, v2) in basis: - Phi_at_qpts = numpy.outer(v1, v2)[None, :, :, None] * phi_at_qpts[:, None, None, :] - nodes.extend(FIM(ref_el, Q, Phi) for Phi in Phi_at_qpts) - - # More internal dofs: tangential-tangential moments against bubbles for each edge - # Note these are evaluated on the edge, but not shared between cells (hence internal). - facet = ref_el.get_facet_element() - Q = create_quadrature(facet, 2*degree-sd) - phi = polynomial_set.ONPolynomialSet(facet, degree-sd) - phi_at_qpts = phi.tabulate(Q.get_points())[(0,) * (sd-1)] - for entity_id in sorted(top[1]): - Q_mapped = FacetQuadratureRule(ref_el, sd-1, entity_id, Q) - t = ref_el.compute_edge_tangent(entity_id) - Phi_at_qpts = numpy.outer(t, t)[None, :, :, None] * phi_at_qpts[:, None, None, :] - nodes.extend(FIM(ref_el, Q_mapped, Phi) for Phi in Phi_at_qpts) - - elif variant == "point": + if variant == "point": # Evaluation at interior points for each component - interior_points = ref_el.make_points(sd, 0, degree) + interior_points = ref_el.make_points(sd, 0, degree+1) nodes.extend(InnerProduct(ref_el, v1, v2, pt) for pt in interior_points for (v1, v2) in basis) - # More internal dofs: tangential-tangential point evaluations. - # Note these are evaluated on the edge, but not shared between cells (hence internal). - for entity_id in sorted(top[1]): - t = ref_el.compute_edge_tangent(entity_id) - pts = ref_el.make_points(sd-1, entity_id, degree) - nodes.extend(InnerProduct(ref_el, t, t, pt) for pt in pts) + elif variant == "integral": + Q = create_quadrature(ref_el, degree + qdegree) + qpts = Q.get_points() + P = polynomial_set.ONPolynomialSet(ref_el, degree-2) + Phis = P.tabulate(qpts)[(0,)*sd] + v = numpy.array(ref_el.vertices) + x = numpy.transpose(ref_el.compute_barycentric_coordinates(qpts)) + for k in sorted(top[1]): + i = (k+1) % (sd+1) + j = (k+2) % (sd+1) + t = v[i] - v[j] + phis = numpy.outer(t, t)[None, :, :, None] * Phis[:, None, None, :] + phis = numpy.multiply(phis, x[i] * x[j], out=phis) + nodes.extend(FIM(ref_el, Q, phi) for phi in phis) + else: raise ValueError(f"Unsupported variant {variant}") @@ -94,12 +81,13 @@ def __init__(self, ref_el, degree, variant): class HuZhang(finite_element.CiarletElement): """The definition of the Hu-Zhang element. """ - def __init__(self, ref_el, degree=3, variant="integral"): + def __init__(self, ref_el, degree=3, variant="point"): if degree < 3: raise ValueError("Hu-Zhang only defined for degree >= 3") if ref_el.shape != TRIANGLE: raise ValueError("Hu-Zhang only defined on triangles") + variant, qdegree = check_format_variant(variant, degree) poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) - dual = HuZhangDual(ref_el, degree, variant) + dual = HuZhangDual(ref_el, degree, variant, qdegree) formdegree = ref_el.get_spatial_dimension() - 1 super().__init__(poly_set, dual, degree, formdegree=formdegree, mapping="double contravariant piola") diff --git a/test/unit/test_hz.py b/test/unit/test_hz.py deleted file mode 100644 index 0a761f40e..000000000 --- a/test/unit/test_hz.py +++ /dev/null @@ -1,136 +0,0 @@ -import pytest -import numpy as np - -from FIAT import ufc_simplex, HuZhang, expansions -from FIAT.quadrature_schemes import create_quadrature - - -def test_dofs(): - line = ufc_simplex(1) - T = ufc_simplex(2) - T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)]) - HZ = HuZhang(T, 3) - - # check Kronecker property at vertices - - bases = [[[1, 0], [0, 0]], [[0, 1], [1, 0]], [[0, 0], [0, 1]]] - - vert_vals = HZ.tabulate(0, T.vertices)[(0, 0)] - for i in range(3): - for j in range(3): - assert np.allclose(vert_vals[3*i+j, :, :, i], bases[j]) - for k in (1, 2): - assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], 0) - - # check edge moments - Qline = create_quadrature(line, 6) - - linebfs = expansions.LineExpansionSet(line) - linevals = linebfs.tabulate(1, Qline.pts) - - # n, n moments - for ed in range(3): - n = T.compute_scaled_normal(ed) - wts = np.asarray(Qline.wts) - nqpline = len(wts) - - vals = HZ.tabulate(0, Qline.pts, (1, ed))[(0, 0)] - nnvals = np.zeros((30, nqpline)) - for i in range(30): - for j in range(len(wts)): - nnvals[i, j] = n @ vals[i, :, :, j] @ n - - nnmoments = np.zeros((30, 2)) - - for bf in range(30): - for k in range(nqpline): - for m in (0, 1): - nnmoments[bf, m] += wts[k] * nnvals[bf, k] * linevals[m, k] - - for bf in range(30): - if bf != HZ.dual.entity_ids[1][ed][0] and bf != HZ.dual.entity_ids[1][ed][2]: - assert np.allclose(nnmoments[bf, :], 0) - - # n, t moments - for ed in range(3): - n = T.compute_scaled_normal(ed) - t = T.compute_edge_tangent(ed) - wts = np.asarray(Qline.wts) - nqpline = len(wts) - - vals = HZ.tabulate(0, Qline.pts, (1, ed))[(0, 0)] - ntvals = np.zeros((30, nqpline)) - for i in range(30): - for j in range(len(wts)): - ntvals[i, j] = n @ vals[i, :, :, j] @ t - - ntmoments = np.zeros((30, 2)) - - for bf in range(30): - for k in range(nqpline): - for m in (0, 1): - ntmoments[bf, m] += wts[k] * ntvals[bf, k] * linevals[m, k] - - for bf in range(30): - if bf != HZ.dual.entity_ids[1][ed][1] and bf != HZ.dual.entity_ids[1][ed][3]: - assert np.allclose(ntmoments[bf, :], 0) - - # check internal dofs - Q = create_quadrature(T, 6) - qpvals = HZ.tabulate(0, Q.pts)[(0, 0)] - const_moms = qpvals @ Q.wts - assert np.allclose(const_moms[:21], 0) - assert np.allclose(const_moms[24:], 0) - assert np.allclose(const_moms[21:24, 0, 0], [1, 0, 0]) - assert np.allclose(const_moms[21:24, 0, 1], [0, 1, 0]) - assert np.allclose(const_moms[21:24, 1, 0], [0, 1, 0]) - assert np.allclose(const_moms[21:24, 1, 1], [0, 0, 1]) - - -def frob(a, b): - return a.ravel() @ b.ravel() - - -@pytest.mark.parametrize("variant", ("integral", "point")) -def test_projection(variant): - T = ufc_simplex(2) - T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)]) - - p = 3 - HZ = HuZhang(T, p, variant) - - Q = create_quadrature(T, 6) - qpts = np.asarray(Q.pts) - qwts = np.asarray(Q.wts) - nqp = len(Q.wts) - - nbf = HZ.space_dimension() - 3 * (p-1) - rhs_vals = np.zeros((1, 2, 2, nqp)) - - bfvals = HZ.tabulate(0, qpts)[(0, 0)][:nbf, :, :, :] - ells = np.multiply(bfvals, qwts) - m = np.tensordot(ells, bfvals, (range(1, ells.ndim),)*2) - - assert np.linalg.cond(m) < 1.e12 - - comps = [(0, 0), (0, 1), (1, 1)] - - # loop over monomials up to degree 2 - for deg in range(3): - for jj in range(deg+1): - ii = deg-jj - for comp in comps: - # set RHS (symmetrically) to be the monomial in - # the proper component. - rhs_vals[...] = 0 - rhs_vals[0][comp] = qpts[:, 0]**ii * qpts[:, 1]**jj - rhs_vals[0][tuple(reversed(comp))] = rhs_vals[0][comp] - - b = np.tensordot(ells, rhs_vals, (range(1, ells.ndim),)*2) - x = np.linalg.solve(m, b) - - sol_at_qpts = np.tensordot(x, bfvals, (0, 0)) - - diff = (sol_at_qpts - rhs_vals)**2 - err = np.linalg.norm(np.tensordot(diff, qwts, (-1, -1))[0], "fro") - assert np.sqrt(err) < 1.e-12 From 66fcb2293a0904a6308ab686dd8f6b336e97bf09 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 12:08:37 +0100 Subject: [PATCH 13/25] Fix HZ integral variant, clean up AW --- FIAT/arnold_winther.py | 222 ++++++++++++++++++----------------------- FIAT/hu_zhang.py | 73 ++++++-------- test/unit/test_awc.py | 4 +- test/unit/test_fiat.py | 1 + 4 files changed, 127 insertions(+), 173 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index 84febd88c..b3aac9539 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -9,154 +9,122 @@ # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT.finite_element import CiarletElement -from FIAT.dual_set import DualSet -from FIAT.polynomial_set import ONSymTensorPolynomialSet, ONPolynomialSet -from FIAT.functional import ( - PointwiseInnerProductEvaluation as InnerProduct, - FrobeniusIntegralMoment as FIM, - IntegralMomentOfTensorDivergence, - IntegralLegendreNormalNormalMoment, - IntegralLegendreNormalTangentialMoment) - -from FIAT.quadrature import make_quadrature +from FIAT import finite_element, dual_set, expansions, polynomial_set +from FIAT.reference_element import TRIANGLE +from FIAT.quadrature_schemes import create_quadrature +from FIAT.functional import (ComponentPointEvaluation, + IntegralMoment, + IntegralMomentOfTensorDivergence, + IntegralLegendreNormalNormalMoment, + IntegralLegendreNormalTangentialMoment) + import numpy -class ArnoldWintherNCDual(DualSet): - def __init__(self, cell, degree=2): - if not degree == 2: - raise ValueError("Nonconforming Arnold-Winther elements are" - "only defined for degree 2.") - dofs = [] - dof_ids = {} - dof_ids[0] = {0: [], 1: [], 2: []} - dof_ids[1] = {0: [], 1: [], 2: []} - dof_ids[2] = {0: []} +class ArnoldWintherNCDual(dual_set.DualSet): + def __init__(self, ref_el, degree=2): + if degree != 2: + raise ValueError("Nonconforming Arnold-Winther elements are only defined for degree 2.") + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + shp = (sd, sd) + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] - dof_cur = 0 # no vertex dofs # proper edge dofs now (not the contraints) - # moments of normal . sigma against constants and linears. - for entity_id in range(3): # a triangle has 3 edges - for order in (0, 1): - dofs += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, degree), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, degree)] - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) - dof_cur += 4 + # edge dofs: bidirectional nn and nt moments against P1. + for entity in sorted(top[1]): + cur = len(nodes) + for order in range(2): + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, degree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, degree)) + entity_ids[1][entity].extend(range(cur, len(nodes))) # internal dofs: constant moments of three unique components - Q = make_quadrature(cell, 2) - - e1 = numpy.array([1.0, 0.0]) # euclidean basis 1 - e2 = numpy.array([0.0, 1.0]) # euclidean basis 2 - basis = [(e1, e1), (e1, e2), (e2, e2)] # basis for symmetric matrices - for (v1, v2) in basis: - v1v2t = numpy.outer(v1, v2) - fatqp = numpy.zeros((2, 2, len(Q.pts))) - for i, y in enumerate(v1v2t): - for j, x in enumerate(y): - for k in range(len(Q.pts)): - fatqp[i, j, k] = x - dofs.append(FIM(cell, Q, fatqp)) - dof_ids[2][0] = list(range(dof_cur, dof_cur + 3)) - dof_cur += 3 + cur = len(nodes) + Q = create_quadrature(ref_el, degree) + phi = numpy.ones(Q.get_weights().shape) + nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + for i in range(sd) for j in range(i, sd)) + entity_ids[2][0].extend(range(cur, len(nodes))) # put the constraint dofs last. - for entity_id in range(3): - dof = IntegralLegendreNormalNormalMoment(cell, entity_id, 2, 6) - dofs.append(dof) - dof_ids[1][entity_id].append(dof_cur) - dof_cur += 1 + for entity in sorted(top[1]): + cur = len(nodes) + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, 2, degree)) + entity_ids[1][entity].append(cur) - super().__init__(dofs, cell, dof_ids) + super().__init__(nodes, ref_el, entity_ids) -class ArnoldWintherNC(CiarletElement): +class ArnoldWintherNC(finite_element.CiarletElement): """The definition of the nonconforming Arnold-Winther element. """ - def __init__(self, cell, degree=2): - assert degree == 2, "Only defined for degree 2" - Ps = ONSymTensorPolynomialSet(cell, degree) - Ls = ArnoldWintherNCDual(cell, degree) + def __init__(self, ref_el, degree=2): + if ref_el.shape != TRIANGLE: + raise ValueError(f"{type(self).__name__} only defined on triangles") + Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + Ls = ArnoldWintherNCDual(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" - - super().__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(Ps, Ls, degree, formdegree, mapping=mapping) -class ArnoldWintherDual(DualSet): - def __init__(self, cell, degree=3): - if not degree == 3: - raise ValueError("Arnold-Winther elements are" - "only defined for degree 3.") - dofs = [] - dof_ids = {} - dof_ids[0] = {0: [], 1: [], 2: []} - dof_ids[1] = {0: [], 1: [], 2: []} - dof_ids[2] = {0: []} - - dof_cur = 0 +class ArnoldWintherDual(dual_set.DualSet): + def __init__(self, ref_el, degree=3): + if degree != 3: + raise ValueError("Arnold-Winther elements are only defined for degree 3.") + top = ref_el.get_topology() + sd = ref_el.get_spatial_dimension() + shp = (sd, sd) + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] # vertex dofs - vs = cell.get_vertices() - e1 = numpy.array([1.0, 0.0]) - e2 = numpy.array([0.0, 1.0]) - basis = [(e1, e1), (e1, e2), (e2, e2)] - - dof_cur = 0 - - for entity_id in range(3): - node = tuple(vs[entity_id]) - for (v1, v2) in basis: - dofs.append(InnerProduct(cell, v1, v2, node)) - dof_ids[0][entity_id] = list(range(dof_cur, dof_cur + 3)) - dof_cur += 3 - - # edge dofs now - # moments of normal . sigma against constants and linears. - for entity_id in range(3): - for order in (0, 1): - dofs.append(IntegralLegendreNormalNormalMoment(cell, entity_id, order, degree)) - dofs.append(IntegralLegendreNormalTangentialMoment(cell, entity_id, order, degree)) - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) - dof_cur += 4 - - # internal dofs: constant moments of three unique components - Q = make_quadrature(cell, 3) - - for (v1, v2) in basis: - v1v2t = numpy.outer(v1, v2) - fatqp = numpy.zeros((2, 2, len(Q.pts))) - for k in range(len(Q.pts)): - fatqp[:, :, k] = v1v2t - dofs.append(FIM(cell, Q, fatqp)) - dof_ids[2][0] = list(range(dof_cur, dof_cur + 3)) - dof_cur += 3 - - # Constraint dofs - - Q = make_quadrature(cell, 5) - - onp = ONPolynomialSet(cell, 2, (2,)) - pts = Q.get_points() - onpvals = onp.tabulate(pts)[0, 0] - - for i in list(range(3, 6)) + list(range(9, 12)): - dofs.append(IntegralMomentOfTensorDivergence(cell, Q, - onpvals[i, :, :])) - - dof_ids[2][0] += list(range(dof_cur, dof_cur+6)) - - super().__init__(dofs, cell, dof_ids) - - -class ArnoldWinther(CiarletElement): + for v in sorted(top[0]): + cur = len(nodes) + pt, = ref_el.make_points(0, v, degree) + nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt) + for i in range(sd) for j in range(i, sd)) + entity_ids[0][v].extend(range(cur, len(nodes))) + + # edge dofs: bidirectional nn and nt moments against P_{k-2} + for entity in sorted(top[1]): + cur = len(nodes) + for order in range(degree-1): + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, degree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, degree)) + entity_ids[1][entity].extend(range(cur, len(nodes))) + + # internal dofs: moments of unique components against P_{k-3} + Q = create_quadrature(ref_el, 2*(degree-1)) + P = polynomial_set.ONPolynomialSet(ref_el, degree-3) + phis = P.tabulate(Q.get_points())[(0,)*sd] + nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + for phi in phis for i in range(sd) for j in range(i, sd)) + + # constraint dofs: moments of divergence against P_{k-1}^\perp P_{k-2} + dimPkm1 = expansions.polynomial_dimension(ref_el, degree-1) + dimPkm2 = expansions.polynomial_dimension(ref_el, degree-2) + P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,)) + PH = P.take([i + j * dimPkm1 for j in range(sd) for i in range(dimPkm2, dimPkm1)]) + phis = PH.tabulate(Q.get_points())[(0,)*sd] + nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) + + entity_ids[2][0].extend(range(cur, len(nodes))) + super().__init__(nodes, ref_el, entity_ids) + + +class ArnoldWinther(finite_element.CiarletElement): """The definition of the conforming Arnold-Winther element. """ - def __init__(self, cell, degree=3): - assert degree == 3, "Only defined for degree 3" - Ps = ONSymTensorPolynomialSet(cell, degree) - Ls = ArnoldWintherDual(cell, degree) + def __init__(self, ref_el, degree=3): + if ref_el.shape != TRIANGLE: + raise ValueError(f"{type(self).__name__} only defined on triangles") + Ps = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) + Ls = ArnoldWintherDual(ref_el, degree) + formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" - super().__init__(Ps, Ls, degree, mapping=mapping) + super().__init__(Ps, Ls, degree, formdegree, mapping=mapping) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 83b486116..58936b941 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -12,13 +12,10 @@ from FIAT.check_format_variant import check_format_variant from FIAT.reference_element import TRIANGLE from FIAT.quadrature_schemes import create_quadrature -from FIAT.functional import (PointwiseInnerProductEvaluation as InnerProduct, - FrobeniusIntegralMoment as FIM, +from FIAT.functional import (ComponentPointEvaluation, + IntegralMoment, IntegralLegendreNormalNormalMoment, - IntegralLegendreNormalTangentialMoment, - ) - -import numpy + IntegralLegendreNormalTangentialMoment) class HuZhangDual(dual_set.DualSet): @@ -27,65 +24,53 @@ def __init__(self, ref_el, degree, variant, qdegree): qdegree = degree top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() + shp = (sd, sd) entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} nodes = [] # vertex dofs - e1 = numpy.array([1.0, 0.0]) - e2 = numpy.array([0.0, 1.0]) - basis = [(e1, e1), (e1, e2), (e2, e2)] - for entity_id in sorted(top[0]): + for v in sorted(top[0]): cur = len(nodes) - pt, = ref_el.make_points(0, entity_id, degree) - nodes.extend(InnerProduct(ref_el, v1, v2, pt) for (v1, v2) in basis) - entity_ids[0][entity_id].extend(range(cur, len(nodes))) + pt, = ref_el.make_points(0, v, degree) + nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt) + for i in range(sd) for j in range(i, sd)) + entity_ids[0][v].extend(range(cur, len(nodes))) - # edge dofs: moments of normal component of sigma against degree p - 2. - for entity_id in sorted(top[1]): + # edge dofs: bidirectional nn and nt moments against P_{k-2}. + for entity in sorted(top[1]): cur = len(nodes) - for k in range(degree-1): - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity_id, k, qdegree)) - nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity_id, k, qdegree)) - entity_ids[1][entity_id].extend(range(cur, len(nodes))) + for order in range(degree-1): + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree)) + entity_ids[1][entity].extend(range(cur, len(nodes))) - # internal dofs + # interior dofs cur = len(nodes) if variant == "point": - # Evaluation at interior points for each component - interior_points = ref_el.make_points(sd, 0, degree+1) - nodes.extend(InnerProduct(ref_el, v1, v2, pt) - for pt in interior_points for (v1, v2) in basis) + # unique components evaluated at interior points + pts = ref_el.make_points(sd, 0, degree+1) + nodes.extend(ComponentPointEvaluation(ref_el, (i, j), shp, pt) + for pt in pts for i in range(sd) for j in range(i, sd)) elif variant == "integral": - Q = create_quadrature(ref_el, degree + qdegree) - qpts = Q.get_points() + # Moments of unique components against a basis for P_{k-2} + Q = create_quadrature(ref_el, qdegree + degree-2) P = polynomial_set.ONPolynomialSet(ref_el, degree-2) - Phis = P.tabulate(qpts)[(0,)*sd] - v = numpy.array(ref_el.vertices) - x = numpy.transpose(ref_el.compute_barycentric_coordinates(qpts)) - for k in sorted(top[1]): - i = (k+1) % (sd+1) - j = (k+2) % (sd+1) - t = v[i] - v[j] - phis = numpy.outer(t, t)[None, :, :, None] * Phis[:, None, None, :] - phis = numpy.multiply(phis, x[i] * x[j], out=phis) - nodes.extend(FIM(ref_el, Q, phi) for phi in phis) - - else: - raise ValueError(f"Unsupported variant {variant}") + phis = P.tabulate(Q.get_points())[(0,)*sd] + nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + for phi in phis for i in range(sd) for j in range(i, sd)) entity_ids[2][0].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) class HuZhang(finite_element.CiarletElement): - """The definition of the Hu-Zhang element. - """ - def __init__(self, ref_el, degree=3, variant="point"): + """The definition of the Hu-Zhang element.""" + def __init__(self, ref_el, degree=3, variant=None): if degree < 3: - raise ValueError("Hu-Zhang only defined for degree >= 3") + raise ValueError(f"{type(self).__name__} only defined for degree >= 3") if ref_el.shape != TRIANGLE: - raise ValueError("Hu-Zhang only defined on triangles") + raise ValueError(f"{type(self).__name__} only defined on triangles") variant, qdegree = check_format_variant(variant, degree) poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) dual = HuZhangDual(ref_el, degree, variant, qdegree) diff --git a/test/unit/test_awc.py b/test/unit/test_awc.py index f40690c48..c5387803a 100644 --- a/test/unit/test_awc.py +++ b/test/unit/test_awc.py @@ -5,7 +5,7 @@ def test_dofs(): line = ufc_simplex(1) T = ufc_simplex(2) - T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)]) + T.vertices = ((0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)) AW = ArnoldWinther(T, 3) # check Kronecker property at vertices @@ -90,7 +90,7 @@ def frob(a, b): def test_projection(): T = ufc_simplex(2) - T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.0), (0.5, 2.1)]) + T.vertices = ((0.0, 0.0), (1.0, 0.0), (0.5, 2.1)) AW = ArnoldWinther(T, 3) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 662212422..389de7ff6 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -310,6 +310,7 @@ def __init__(self, a, b): "Morley(T)", "BernardiRaugel(T)", "BernardiRaugel(S)", + "ArnoldWintherNC(T, 2)", "ArnoldWinther(T, 3)", "HuZhang(T, 3)", "HuZhang(T, 4)", From c800a5d37ef823f288501a57d2b5e1dd56ca3b7f Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 14:36:24 +0100 Subject: [PATCH 14/25] cache quadrature rules --- FIAT/__init__.py | 12 ++--- FIAT/functional.py | 90 ++++++++++++++------------------------ FIAT/quadrature_schemes.py | 15 +++++-- 3 files changed, 50 insertions(+), 67 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 8b522ee2d..c56daf965 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -4,6 +4,12 @@ import pkg_resources +# Important functionality +from FIAT.quadrature import make_quadrature # noqa: F401 +from FIAT.quadrature_schemes import create_quadrature # noqa: F401 +from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401 +from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 + # Import finite element classes from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 from FIAT.argyris import Argyris @@ -59,12 +65,6 @@ from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401 from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401 -# Important functionality -from FIAT.quadrature import make_quadrature # noqa: F401 -from FIAT.quadrature_schemes import create_quadrature # noqa: F401 -from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401 -from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 - __version__ = pkg_resources.get_distribution("fenics-fiat").version # List of supported elements and mapping to element classes diff --git a/FIAT/functional.py b/FIAT/functional.py index 3cc5c88de..865598d21 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -16,9 +16,7 @@ import numpy import sympy -from FIAT import polynomial_set, jacobi -from FIAT.quadrature import GaussLegendreQuadratureLineRule -from FIAT.reference_element import UFCInterval as interval +from FIAT import polynomial_set, jacobi, quadrature_schemes def index_iterator(shp): @@ -338,23 +336,39 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts): {}, dpt_dict, "IntegralMomentOfNormalDerivative") -class IntegralLegendreDirectionalMoment(Functional): +class FrobeniusIntegralMoment(IntegralMoment): + + def __init__(self, ref_el, Q, f_at_qpts): + # f_at_qpts is (some shape) x num_qpts + shp = tuple(f_at_qpts.shape[:-1]) + if len(Q.pts) != f_at_qpts.shape[-1]: + raise Exception("Mismatch in number of quadrature points and values") + + self.Q = Q + self.comp = slice(None, None) + self.f_at_qpts = f_at_qpts + qpts, qwts = Q.get_points(), Q.get_weights() + weights = numpy.transpose(numpy.multiply(f_at_qpts, qwts), (-1,) + tuple(range(len(shp)))) + alphas = list(index_iterator(shp)) + + pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} + Functional.__init__(self, ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment") + + +class IntegralLegendreDirectionalMoment(FrobeniusIntegralMoment): """Moment of v.s against a Legendre polynomial over an edge""" def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): - sd = cell.get_spatial_dimension() - assert sd == 2 - shp = (sd,) - quadpoints = comp_deg + 1 - Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) - x = 2*Q.get_points()[:, 0]-1 - f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x) - transform = cell.get_entity_transform(sd-1, entity) - points = transform(Q.get_points()) - weights = numpy.multiply(f_at_qpts, Q.get_weights()) - pt_dict = {tuple(pt): [(wt*s[i], (i,)) for i in range(sd)] - for pt, wt in zip(points, weights)} + assert cell.get_spatial_dimension() == 2 + entity = (1, entity) - super().__init__(cell, shp, pt_dict, {}, nm) + Q = quadrature_schemes.create_quadrature(cell, 2*comp_deg, entity=entity) + x = cell.compute_barycentric_coordinates(Q.get_points(), entity=entity) + + f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x[:, 1] - x[:, 0]) + f_at_qpts /= Q.jacobian_determinant() + + f_at_qpts = numpy.multiply(s[..., None], f_at_qpts) + super().__init__(cell, Q, f_at_qpts) class IntegralLegendreNormalMoment(IntegralLegendreDirectionalMoment): @@ -373,32 +387,13 @@ def __init__(self, cell, entity, mom_deg, comp_deg): "IntegralLegendreTangentialMoment") -class IntegralLegendreBidirectionalMoment(Functional): +class IntegralLegendreBidirectionalMoment(IntegralLegendreDirectionalMoment): """Moment of dot(s1, dot(tau, s2)) against Legendre on entity, multiplied by the size of the reference facet""" def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""): # mom_deg is degree of moment, comp_deg is the total degree of # polynomial you might need to integrate (or something like that) - sd = cell.get_spatial_dimension() - s1s2T = numpy.outer(s1, s2) - shp = s1s2T.shape - quadpoints = comp_deg + 1 - Q = GaussLegendreQuadratureLineRule(interval(), quadpoints) - - # The volume squared gets the Jacobian mapping from line interval - # and the edge length into the functional. - x = 2*Q.get_points()[:, 0]-1 - f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x) * numpy.abs(cell.volume_of_subcomplex(1, entity))**2 - - # Map the quadrature points - transform = cell.get_entity_transform(sd-1, entity) - points = transform(Q.get_points()) - weights = numpy.multiply(f_at_qpts, Q.get_weights()) - - pt_dict = {tuple(pt): [(wt * s1s2T[idx], idx) for idx in index_iterator(shp)] - for pt, wt in zip(points, weights)} - - super().__init__(cell, shp, pt_dict, {}, nm) + super().__init__(cell, s1s2T, entity, mom_deg, comp_deg, nm=nm) class IntegralLegendreNormalNormalMoment(IntegralLegendreBidirectionalMoment): @@ -470,25 +465,6 @@ def __init__(self, ref_el, Q, f_at_qpts): super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence") -class FrobeniusIntegralMoment(IntegralMoment): - - def __init__(self, ref_el, Q, f_at_qpts): - # f_at_qpts is (some shape) x num_qpts - shp = tuple(f_at_qpts.shape[:-1]) - if len(Q.pts) != f_at_qpts.shape[-1]: - raise Exception("Mismatch in number of quadrature points and values") - - self.Q = Q - self.comp = slice(None, None) - self.f_at_qpts = f_at_qpts - qpts, qwts = Q.get_points(), Q.get_weights() - weights = numpy.transpose(numpy.multiply(f_at_qpts, qwts), (-1,) + tuple(range(len(shp)))) - alphas = list(index_iterator(shp)) - - pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} - Functional.__init__(self, ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment") - - class PointNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a point on a facet of codimension 1.""" diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index 578c9800a..ca58d67f6 100644 --- a/FIAT/quadrature_schemes.py +++ b/FIAT/quadrature_schemes.py @@ -29,17 +29,18 @@ # NumPy import numpy +import functools -from FIAT.macro import MacroQuadratureRule -from FIAT.quadrature import (QuadratureRule, make_quadrature, +from FIAT.quadrature import (QuadratureRule, FacetQuadratureRule, make_quadrature, make_tensor_product_quadrature, map_quadrature) -# FIAT from FIAT.reference_element import (HEXAHEDRON, QUADRILATERAL, TENSORPRODUCT, TETRAHEDRON, TRIANGLE, UFCTetrahedron, UFCTriangle, symmetric_simplex) +from FIAT.macro import MacroQuadratureRule -def create_quadrature(ref_el, degree, scheme="default"): +@functools.lru_cache +def create_quadrature(ref_el, degree, scheme="default", entity=None): """ Generate quadrature rule for given reference element that will integrate an polynomial of order 'degree' exactly. @@ -53,6 +54,12 @@ def create_quadrature(ref_el, degree, scheme="default"): :arg degree: The degree of polynomial that the rule should integrate exactly. """ + if entity is not None: + dim, entity_id = entity + sub_el = ref_el.construct_subelement(dim) + Q_ref = create_quadrature(sub_el, degree, scheme=scheme) + return FacetQuadratureRule(ref_el, dim, entity_id, Q_ref) + if ref_el.is_macrocell(): dimension = ref_el.get_dimension() sub_el = ref_el.construct_subelement(dimension) From 8bf370890abd64fa8f5b5a6a3ed9434a1abafa45 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 16:32:09 +0100 Subject: [PATCH 15/25] docs --- FIAT/functional.py | 12 ++++++------ FIAT/quadrature_schemes.py | 23 ++++++++++++++--------- 2 files changed, 20 insertions(+), 15 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index 865598d21..066de3786 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -191,7 +191,7 @@ def __init__(self, ref_el, comp, shp, x): if any(i < 0 or i >= n for i, n in zip(comp, shp)): raise ValueError("Illegal component") self.comp = comp - pt_dict = {x: [(1.0, comp)]} + pt_dict = {tuple(x): [(1.0, comp)]} super().__init__(ref_el, shp, pt_dict, {}, "ComponentPointEval") def tostr(self): @@ -338,7 +338,7 @@ def __init__(self, ref_el, facet_no, Q, f_at_qpts): class FrobeniusIntegralMoment(IntegralMoment): - def __init__(self, ref_el, Q, f_at_qpts): + def __init__(self, ref_el, Q, f_at_qpts, nm=None): # f_at_qpts is (some shape) x num_qpts shp = tuple(f_at_qpts.shape[:-1]) if len(Q.pts) != f_at_qpts.shape[-1]: @@ -352,12 +352,14 @@ def __init__(self, ref_el, Q, f_at_qpts): alphas = list(index_iterator(shp)) pt_dict = {tuple(pt): [(wt[alpha], alpha) for alpha in alphas] for pt, wt in zip(qpts, weights)} - Functional.__init__(self, ref_el, shp, pt_dict, {}, "FrobeniusIntegralMoment") + Functional.__init__(self, ref_el, shp, pt_dict, {}, nm or "FrobeniusIntegralMoment") class IntegralLegendreDirectionalMoment(FrobeniusIntegralMoment): """Moment of v.s against a Legendre polynomial over an edge""" def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): + # mom_deg is degree of moment, comp_deg is the total degree of + # polynomial you might need to integrate (or something like that) assert cell.get_spatial_dimension() == 2 entity = (1, entity) @@ -368,7 +370,7 @@ def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): f_at_qpts /= Q.jacobian_determinant() f_at_qpts = numpy.multiply(s[..., None], f_at_qpts) - super().__init__(cell, Q, f_at_qpts) + super().__init__(cell, Q, f_at_qpts, nm=nm) class IntegralLegendreNormalMoment(IntegralLegendreDirectionalMoment): @@ -390,8 +392,6 @@ def __init__(self, cell, entity, mom_deg, comp_deg): class IntegralLegendreBidirectionalMoment(IntegralLegendreDirectionalMoment): """Moment of dot(s1, dot(tau, s2)) against Legendre on entity, multiplied by the size of the reference facet""" def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""): - # mom_deg is degree of moment, comp_deg is the total degree of - # polynomial you might need to integrate (or something like that) s1s2T = numpy.outer(s1, s2) super().__init__(cell, s1s2T, entity, mom_deg, comp_deg, nm=nm) diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index ca58d67f6..1dff01213 100644 --- a/FIAT/quadrature_schemes.py +++ b/FIAT/quadrature_schemes.py @@ -27,7 +27,6 @@ # First added: 2011-04-19 # Last changed: 2011-04-19 -# NumPy import numpy import functools @@ -42,17 +41,23 @@ @functools.lru_cache def create_quadrature(ref_el, degree, scheme="default", entity=None): """ - Generate quadrature rule for given reference element - that will integrate an polynomial of order 'degree' exactly. + Generate quadrature rule for given reference element that will integrate an + polynomial of order 'degree' exactly. - For low-degree (<=6) polynomials on triangles and tetrahedra, this - uses hard-coded rules, otherwise it falls back to a collapsed - Gauss scheme on simplices. On tensor-product cells, it is a - tensor-product quadrature rule of the subcells. + For low-degree polynomials on triangles (<=50) and tetrahedra (<=15), this uses + hard-coded rules, otherwise it falls back to a collapsed Gauss scheme on + simplices. On tensor-product cells, it is a tensor-product quadrature rule + of the subcells. :arg ref_el: The FIAT cell to create the quadrature for. - :arg degree: The degree of polynomial that the rule should - integrate exactly. + :arg degree: The degree of polynomial that the rule should integrate exactly. + :kwarg scheme: The quadrature scheme, can be choosen from ["default", "canonical", "KMV"] + "default" -> optimized Xiao-Gimbutas scheme for low degree and + collapsed Gauss scheme for higher degree, + "canonical" -> collapsed Gauss scheme, + "KMV" -> spectral lumped scheme for low degree (<=5 on triangles, <=3 on tetrahedra). + :kwarg entity: A tuple of entity dimension and entity id specifying the + integration domain. If not provided, the domain is the entire cell. """ if entity is not None: dim, entity_id = entity From 417466eb87c0aeb7a50f6d58d2fd801acce1d2e5 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 19:43:09 +0100 Subject: [PATCH 16/25] Fix quadrature degrees, cleanup MTW --- FIAT/arnold_winther.py | 17 +++++++------ FIAT/functional.py | 38 ++++++----------------------- FIAT/hu_zhang.py | 21 +++++++++++----- FIAT/mardal_tai_winther.py | 50 ++++++++++++++++---------------------- FIAT/reference_element.py | 2 +- test/unit/test_fiat.py | 2 ++ 6 files changed, 57 insertions(+), 73 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index b3aac9539..aa2523e9a 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -35,11 +35,12 @@ def __init__(self, ref_el, degree=2): # no vertex dofs # proper edge dofs now (not the contraints) # edge dofs: bidirectional nn and nt moments against P1. + qdegree = degree + 2 for entity in sorted(top[1]): cur = len(nodes) for order in range(2): - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, degree)) - nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, degree)) + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree)) entity_ids[1][entity].extend(range(cur, len(nodes))) # internal dofs: constant moments of three unique components @@ -53,7 +54,7 @@ def __init__(self, ref_el, degree=2): # put the constraint dofs last. for entity in sorted(top[1]): cur = len(nodes) - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, 2, degree)) + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, 2, qdegree)) entity_ids[1][entity].append(cur) super().__init__(nodes, ref_el, entity_ids) @@ -91,11 +92,13 @@ def __init__(self, ref_el, degree=3): entity_ids[0][v].extend(range(cur, len(nodes))) # edge dofs: bidirectional nn and nt moments against P_{k-2} + max_order = degree - 2 + qdegree = degree + max_order for entity in sorted(top[1]): cur = len(nodes) - for order in range(degree-1): - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, degree)) - nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, degree)) + for order in range(max_order+1): + nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree)) + nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree)) entity_ids[1][entity].extend(range(cur, len(nodes))) # internal dofs: moments of unique components against P_{k-3} @@ -105,7 +108,7 @@ def __init__(self, ref_el, degree=3): nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) for phi in phis for i in range(sd) for j in range(i, sd)) - # constraint dofs: moments of divergence against P_{k-1}^\perp P_{k-2} + # constraint dofs: moments of divergence against P_{k-1} \ P_{k-2} dimPkm1 = expansions.polynomial_dimension(ref_el, degree-1) dimPkm2 = expansions.polynomial_dimension(ref_el, degree-2) P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,)) diff --git a/FIAT/functional.py b/FIAT/functional.py index 066de3786..629da0ca5 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -294,10 +294,10 @@ class IntegralMoment(Functional): def __init__(self, ref_el, Q, f_at_qpts, comp=tuple(), shp=tuple()): self.Q = Q self.f_at_qpts = f_at_qpts - qpts, qwts = Q.get_points(), Q.get_weights() self.comp = comp - weights = numpy.multiply(f_at_qpts, qwts) - pt_dict = {tuple(pt): [(wt, comp)] for pt, wt in zip(qpts, weights)} + points = Q.get_points() + weights = numpy.multiply(f_at_qpts, Q.get_weights()) + pt_dict = {tuple(pt): [(wt, comp)] for pt, wt in zip(points, weights)} super().__init__(ref_el, shp, pt_dict, {}, "IntegralMoment") def __call__(self, fn): @@ -357,13 +357,13 @@ def __init__(self, ref_el, Q, f_at_qpts, nm=None): class IntegralLegendreDirectionalMoment(FrobeniusIntegralMoment): """Moment of v.s against a Legendre polynomial over an edge""" - def __init__(self, cell, s, entity, mom_deg, comp_deg, nm=""): - # mom_deg is degree of moment, comp_deg is the total degree of + def __init__(self, cell, s, entity, mom_deg, quad_deg, nm=""): + # mom_deg is degree of moment, quad_deg is the total degree of # polynomial you might need to integrate (or something like that) assert cell.get_spatial_dimension() == 2 entity = (1, entity) - Q = quadrature_schemes.create_quadrature(cell, 2*comp_deg, entity=entity) + Q = quadrature_schemes.create_quadrature(cell, quad_deg, entity=entity) x = cell.compute_barycentric_coordinates(Q.get_points(), entity=entity) f_at_qpts = jacobi.eval_jacobi(0, 0, mom_deg, x[:, 1] - x[:, 0]) @@ -565,39 +565,17 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): "IntegralMomentOfFaceTangentEvaluation") -class MonkIntegralMoment(Functional): - r""" - face nodes are \int_F v\cdot p dA where p \in P_{q-2}(f)^3 with p \cdot n = 0 - (cmp. Peter Monk - Finite Element Methods for Maxwell's equations p. 129) - Note that we don't scale by the area of the facet - - :arg ref_el: reference element for which F is a codim-1 entity - :arg Q: quadrature rule on the face - :arg P_at_qpts: polynomials evaluated at quad points - :arg facet: which facet. - """ - - def __init__(self, ref_el, Q, P_at_qpts, facet): - sd = ref_el.get_spatial_dimension() - transform = ref_el.get_entity_transform(sd-1, facet) - pts = transform(Q.get_points()) - weights = Q.get_weights() * P_at_qpts - pt_dict = {tuple(pt): [(wt[i], (i, )) for i in range(sd)] - for pt, wt in zip(pts, weights)} - super().__init__(ref_el, (sd, ), pt_dict, {}, "MonkIntegralMoment") - - class PointScaledNormalEvaluation(Functional): """Implements the evaluation of the normal component of a vector at a point on a facet of codimension 1, where the normal is scaled by the volume of that facet.""" def __init__(self, ref_el, facet_no, pt): - self.n = ref_el.compute_scaled_normal(facet_no) + n = ref_el.compute_scaled_normal(facet_no) sd = ref_el.get_spatial_dimension() shp = (sd,) - pt_dict = {pt: [(self.n[i], (i,)) for i in range(sd)]} + pt_dict = {pt: [(n[i], (i,)) for i in range(sd)]} super().__init__(ref_el, shp, pt_dict, {}, "PointScaledNormalEval") def tostr(self): diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 58936b941..ce7364e45 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -13,6 +13,7 @@ from FIAT.reference_element import TRIANGLE from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (ComponentPointEvaluation, + PointwiseInnerProductEvaluation, IntegralMoment, IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment) @@ -20,8 +21,6 @@ class HuZhangDual(dual_set.DualSet): def __init__(self, ref_el, degree, variant, qdegree): - if qdegree is None: - qdegree = degree top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() shp = (sd, sd) @@ -36,12 +35,22 @@ def __init__(self, ref_el, degree, variant, qdegree): for i in range(sd) for j in range(i, sd)) entity_ids[0][v].extend(range(cur, len(nodes))) - # edge dofs: bidirectional nn and nt moments against P_{k-2}. + # edge dofs for entity in sorted(top[1]): cur = len(nodes) - for order in range(degree-1): - nodes.append(IntegralLegendreNormalNormalMoment(ref_el, entity, order, qdegree)) - nodes.append(IntegralLegendreNormalTangentialMoment(ref_el, entity, order, qdegree)) + if variant == "point": + # nn and nt components evaluated at edge points + n = ref_el.compute_scaled_normal(entity) + t = ref_el.compute_edge_tangent(entity) + pts = ref_el.make_points(1, entity, degree) + nodes.extend(PointwiseInnerProductEvaluation(ref_el, n, s, pt) + for pt in pts for s in (n, t)) + + elif variant == "integral": + # bidirectional nn and nt moments against P_{k-2} + moments = (IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment) + nodes.extend(mu(ref_el, entity, order, qdegree + degree-2) + for order in range(degree-1) for mu in moments) entity_ids[1][entity].extend(range(cur, len(nodes))) # interior dofs diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 0dbc76262..354a71601 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -8,34 +8,28 @@ # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT.finite_element import CiarletElement -from FIAT.dual_set import DualSet -from FIAT.polynomial_set import ONPolynomialSet +from FIAT import dual_set, expansions, finite_element, polynomial_set from FIAT.functional import (IntegralMomentOfNormalEvaluation, IntegralMomentOfTangentialEvaluation, IntegralLegendreNormalMoment, IntegralMomentOfDivergence) -from FIAT.quadrature import make_quadrature +from FIAT.quadrature_schemes import create_quadrature def DivergenceDubinerMoments(cell, start_deg, stop_deg, comp_deg): - onp = ONPolynomialSet(cell, stop_deg) - Q = make_quadrature(cell, comp_deg) + sd = cell.get_spatial_dimension() + P = polynomial_set.ONPolynomialSet(cell, stop_deg) + Q = create_quadrature(cell, comp_deg + stop_deg) - pts = Q.get_points() - onp = onp.tabulate(pts, 0)[0, 0] + dim0 = expansions.polynomial_dimension(cell, start_deg-1) + dim1 = expansions.polynomial_dimension(cell, stop_deg) + indices = list(range(dim0, dim1)) + phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd] + return [IntegralMomentOfDivergence(cell, Q, phi) for phi in phis] - ells = [] - for ii in range((start_deg)*(start_deg+1)//2, - (stop_deg+1)*(stop_deg+2)//2): - ells.append(IntegralMomentOfDivergence(cell, Q, onp[ii, :])) - - return ells - - -class MardalTaiWintherDual(DualSet): +class MardalTaiWintherDual(dual_set.DualSet): """Degrees of freedom for Mardal-Tai-Winther elements.""" def __init__(self, cell, degree): dim = cell.get_spatial_dimension() @@ -93,15 +87,13 @@ def _generate_edge_dofs(cell, degree): facet = cell.get_facet_element() # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1} # degree is q - 1 - Q = make_quadrature(facet, 6) - Pq = ONPolynomialSet(facet, 1) - Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([0]*(sd - 1))] + Q = create_quadrature(facet, degree+1) + Pq = polynomial_set.ONPolynomialSet(facet, 1) + phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)] for f in range(3): - phi0 = Pq_at_qpts[0, :] - dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phi0, f)) - dofs.append(IntegralMomentOfTangentialEvaluation(cell, Q, phi0, f)) - phi1 = Pq_at_qpts[1, :] - dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phi1, f)) + dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phis[0], f)) + dofs.append(IntegralMomentOfTangentialEvaluation(cell, Q, phis[0], f)) + dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phis[1], f)) num_new_dofs = 3 dof_ids[f] = list(range(offset, offset + num_new_dofs)) @@ -129,8 +121,8 @@ def _generate_constraint_dofs(cell, degree, offset): edge_dof_ids = {} for entity_id in range(3): - dofs += [IntegralLegendreNormalMoment(cell, entity_id, 2, 6), - IntegralLegendreNormalMoment(cell, entity_id, 3, 6)] + dofs.append(IntegralLegendreNormalMoment(cell, entity_id, 2, degree+3)) + dofs.append(IntegralLegendreNormalMoment(cell, entity_id, 3, degree+3)) edge_dof_ids[entity_id] = [offset, offset+1] offset += 2 @@ -142,14 +134,14 @@ def _generate_constraint_dofs(cell, degree, offset): return (dofs, edge_dof_ids, cell_dof_ids) -class MardalTaiWinther(CiarletElement): +class MardalTaiWinther(finite_element.CiarletElement): """The definition of the Mardal-Tai-Winther element. """ def __init__(self, cell, degree=3): assert degree == 3, "Only defined for degree 3" assert cell.get_spatial_dimension() == 2, "Only defined for dimension 2" # polynomial space - Ps = ONPolynomialSet(cell, degree, (2,)) + Ps = polynomial_set.ONPolynomialSet(cell, degree, (2,)) # degrees of freedom Ls = MardalTaiWintherDual(cell, degree) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index f09751578..742df82e3 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -570,7 +570,7 @@ def compute_barycentric_coordinates(self, points, entity=None, rescale=False): def compute_bubble(self, points, entity=None): """Returns the lowest-order bubble on an entity evaluated at the given - points on the entity.""" + points on the cell.""" return numpy.prod(self.compute_barycentric_coordinates(points, entity), axis=1) def distance_to_point_l1(self, points, entity=None, rescale=False): diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 389de7ff6..cdf540b2a 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -42,6 +42,7 @@ from FIAT.tensor_product import TensorProductElement # noqa: F401 from FIAT.tensor_product import FlattenedDimensions # noqa: F401 from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 +from FIAT.mardal_tai_winther import MardalTaiWinther # noqa: F401 from FIAT.arnold_winther import ArnoldWinther, ArnoldWintherNC # noqa: F401 from FIAT.hu_zhang import HuZhang # noqa: F401 from FIAT.bernardi_raugel import BernardiRaugel # noqa: F401 @@ -310,6 +311,7 @@ def __init__(self, a, b): "Morley(T)", "BernardiRaugel(T)", "BernardiRaugel(S)", + "MardalTaiWinther(T, 3)", "ArnoldWintherNC(T, 2)", "ArnoldWinther(T, 3)", "HuZhang(T, 3)", From 7397100708e7f80ae7d292e8dafb5916aaf812b1 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sat, 26 Oct 2024 19:49:04 +0100 Subject: [PATCH 17/25] delete unused class --- FIAT/functional.py | 48 ---------------------------------------------- 1 file changed, 48 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index 629da0ca5..30e64252a 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -625,30 +625,6 @@ def __init__(self, ref_el, v, w, pt): super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval") -class TensorBidirectionalMomentInnerProductEvaluation(Functional): - r""" - This is a functional on symmetric 2-tensor fields. Let u be such a - field, f a function tabulated at points, and v,w be vectors. This implements the evaluation - \int v^T u(x) w f(x). - - Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed - from the Frobenius inner product of u with wv^T. This gives the - correct weights. - """ - - def __init__(self, ref_el, v, w, Q, f_at_qpts, comp_deg): - wvT = numpy.outer(w, v) - shp = wvT.shp - - points = Q.get_points() - weights = numpy.multiply(f_at_qpts, Q.get_weights()) - - pt_dict = {tuple(pt): [(wt * wvT[idx], idx) for idx in index_iterator(shp)] - for pt, wt in zip(points, weights)} - - super().__init__(ref_el, shp, pt_dict, {}, "TensorBidirectionalMomentInnerProductEvaluation") - - class IntegralMomentOfNormalEvaluation(Functional): r""" \int_F v\cdot n p ds @@ -692,27 +668,3 @@ def __init__(self, ref_el, Q, P_at_qpts, facet): pt_dict = {tuple(pt): [(wt*t[i], (i, )) for i in range(sd)] for pt, wt in zip(points, weights)} super().__init__(ref_el, (sd, ), pt_dict, {}, "IntegralMomentOfScaledTangentialEvaluation") - - -class IntegralMomentOfNormalNormalEvaluation(Functional): - r""" - \int_F (n^T tau n) p ds - p \in Polynomials - :arg ref_el: reference element for which F is a codim-1 entity - :arg Q: quadrature rule on the face - :arg P_at_qpts: polynomials evaluated at quad points - :arg facet: which facet. - """ - def __init__(self, ref_el, Q, P_at_qpts, facet): - # scaling on the normal is ok because edge length then weights - # the reference element quadrature appropriately - n = ref_el.compute_scaled_normal(facet) - nnT = numpy.outer(n, n)/numpy.linalg.norm(n) - shp = nnT.shape - sd = ref_el.get_spatial_dimension() - transform = ref_el.get_entity_transform(sd - 1, facet) - points = transform(Q.get_points()) - weights = numpy.multiply(P_at_qpts, Q.get_weights()) - pt_dict = {tuple(pt): [(wt*nnT[idx], idx) for idx in index_iterator(shp)] - for pt, wt in zip(points, weights)} - super().__init__(ref_el, shp, pt_dict, {}, "IntegralMomentOfNormalNormalEvaluation") From 6830ed93c5b903d01c18d89b6c8ace3f85cbf026 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 27 Oct 2024 08:36:00 +0000 Subject: [PATCH 18/25] More cleanup --- FIAT/__init__.py | 19 +++++++------------ FIAT/functional.py | 16 ++++++++++++++++ FIAT/quadrature_schemes.py | 16 +++++++++------- 3 files changed, 32 insertions(+), 19 deletions(-) diff --git a/FIAT/__init__.py b/FIAT/__init__.py index c56daf965..47ba99093 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -2,16 +2,17 @@ evaluating arbitrary order Lagrange and many other elements. Simplices in one, two, and three dimensions are supported.""" -import pkg_resources - # Important functionality +from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401 from FIAT.quadrature import make_quadrature # noqa: F401 from FIAT.quadrature_schemes import create_quadrature # noqa: F401 -from FIAT.reference_element import ufc_cell, ufc_simplex # noqa: F401 from FIAT.hdivcurl import Hdiv, Hcurl # noqa: F401 +from FIAT.mixed import MixedElement # noqa: F401 +from FIAT.restricted import RestrictedElement # noqa: F401 +from FIAT.quadrature_element import QuadratureElement # noqa: F401 +from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 # Import finite element classes -from FIAT.finite_element import FiniteElement, CiarletElement # noqa: F401 from FIAT.argyris import Argyris from FIAT.bernardi_raugel import BernardiRaugel from FIAT.bernstein import Bernstein @@ -23,8 +24,7 @@ from FIAT.christiansen_hu import ChristiansenHu from FIAT.johnson_mercier import JohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini -from FIAT.Sminus import TrimmedSerendipityEdge # noqa: F401 -from FIAT.Sminus import TrimmedSerendipityFace # noqa: F401 +from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace # noqa: F401 from FIAT.SminusDiv import TrimmedSerendipityDiv # noqa: F401 from FIAT.SminusCurl import TrimmedSerendipityCurl # noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini @@ -59,14 +59,9 @@ from FIAT.nodal_enriched import NodalEnrichedElement from FIAT.discontinuous import DiscontinuousElement from FIAT.hdiv_trace import HDivTrace -from FIAT.mixed import MixedElement # noqa: F401 -from FIAT.restricted import RestrictedElement # noqa: F401 -from FIAT.quadrature_element import QuadratureElement # noqa: F401 -from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen # noqa: F401 +from FIAT.kong_mulder_veldhuizen import KongMulderVeldhuizen from FIAT.fdm_element import FDMLagrange, FDMDiscontinuousLagrange, FDMQuadrature, FDMBrokenH1, FDMBrokenL2, FDMHermite # noqa: F401 -__version__ = pkg_resources.get_distribution("fenics-fiat").version - # List of supported elements and mapping to element classes supported_elements = {"Argyris": Argyris, "Bell": Bell, diff --git a/FIAT/functional.py b/FIAT/functional.py index 30e64252a..bd0f86889 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -625,6 +625,22 @@ def __init__(self, ref_el, v, w, pt): super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval") +class TensorBidirectionalMomentInnerProductEvaluation(FrobeniusIntegralMoment): + r""" + This is a functional on symmetric 2-tensor fields. Let u be such a + field, f a function tabulated at points, and v,w be vectors. This implements the evaluation + \int v^T u(x) w f(x). + Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed + from the Frobenius inner product of u with wv^T. This gives the + correct weights. + """ + + def __init__(self, ref_el, v, w, Q, f_at_qpts): + wvT = numpy.outer(w, v) + F_at_qpts = numpy.multiply(wvT[..., None], f_at_qpts) + super().__init__(ref_el, Q, F_at_qpts, "TensorBidirectionalMomentInnerProductEvaluation") + + class IntegralMomentOfNormalEvaluation(Functional): r""" \int_F v\cdot n p ds diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index 1dff01213..6449748cb 100644 --- a/FIAT/quadrature_schemes.py +++ b/FIAT/quadrature_schemes.py @@ -15,6 +15,11 @@ Keast, P. Moderate-degree tetrahedral quadrature formulas, Computer Methods in Applied Mechanics and Engineering 55(3):339-348, 1986. http://dx.doi.org/10.1016/0045-7825(86)90059-9 + + Xiao-Gimbutas rules for simplices: + Xiao, H., and Gimbutas, Z. A numerical algorithm for the construction of + efficient quadrature rules in two and higher dimensions, Computers & + mathematics with applications 59(2): 663-676, 2010. """ # Copyright (C) 2011 Garth N. Wells @@ -28,7 +33,6 @@ # Last changed: 2011-04-19 import numpy -import functools from FIAT.quadrature import (QuadratureRule, FacetQuadratureRule, make_quadrature, make_tensor_product_quadrature, map_quadrature) @@ -38,10 +42,9 @@ from FIAT.macro import MacroQuadratureRule -@functools.lru_cache def create_quadrature(ref_el, degree, scheme="default", entity=None): """ - Generate quadrature rule for given reference element that will integrate an + Generate quadrature rule for given reference element that will integrate a polynomial of order 'degree' exactly. For low-degree polynomials on triangles (<=50) and tetrahedra (<=15), this uses @@ -52,8 +55,7 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None): :arg ref_el: The FIAT cell to create the quadrature for. :arg degree: The degree of polynomial that the rule should integrate exactly. :kwarg scheme: The quadrature scheme, can be choosen from ["default", "canonical", "KMV"] - "default" -> optimized Xiao-Gimbutas scheme for low degree and - collapsed Gauss scheme for higher degree, + "default" -> hard-coded scheme for low degree and collapsed Gauss scheme for high degree, "canonical" -> collapsed Gauss scheme, "KMV" -> spectral lumped scheme for low degree (<=5 on triangles, <=3 on tetrahedra). :kwarg entity: A tuple of entity dimension and entity id specifying the @@ -66,8 +68,8 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None): return FacetQuadratureRule(ref_el, dim, entity_id, Q_ref) if ref_el.is_macrocell(): - dimension = ref_el.get_dimension() - sub_el = ref_el.construct_subelement(dimension) + dim = ref_el.get_dimension() + sub_el = ref_el.construct_subelement(dim) Q_ref = create_quadrature(sub_el, degree, scheme=scheme) return MacroQuadratureRule(ref_el, Q_ref) From 2075f9b9da8ae90d38d7ff8d08547a001704d699 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 27 Oct 2024 08:57:17 +0000 Subject: [PATCH 19/25] fix scale --- FIAT/functional.py | 8 ++++---- FIAT/quadrature_schemes.py | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/FIAT/functional.py b/FIAT/functional.py index bd0f86889..817701187 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -399,7 +399,7 @@ def __init__(self, cell, s1, s2, entity, mom_deg, comp_deg, nm=""): class IntegralLegendreNormalNormalMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(n, dot(tau, n)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): - n = cell.compute_normal(entity) + n = cell.compute_scaled_normal(entity) super().__init__(cell, n, n, entity, mom_deg, comp_deg, "IntegralNormalNormalLegendreMoment") @@ -407,8 +407,8 @@ def __init__(self, cell, entity, mom_deg, comp_deg): class IntegralLegendreNormalTangentialMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(n, dot(tau, t)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): - n = cell.compute_normal(entity) - t = cell.compute_normalized_edge_tangent(entity) + n = cell.compute_scaled_normal(entity) + t = cell.compute_edge_tangent(entity) super().__init__(cell, n, t, entity, mom_deg, comp_deg, "IntegralNormalTangentialLegendreMoment") @@ -416,7 +416,7 @@ def __init__(self, cell, entity, mom_deg, comp_deg): class IntegralLegendreTangentialTangentialMoment(IntegralLegendreBidirectionalMoment): """Moment of dot(t, dot(tau, t)) against Legendre on entity.""" def __init__(self, cell, entity, mom_deg, comp_deg): - t = cell.compute_normalized_edge_tangent(entity) + t = cell.compute_edge_tangent(entity) super().__init__(cell, t, t, entity, mom_deg, comp_deg, "IntegralTangentialTangentialLegendreMoment") diff --git a/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index 6449748cb..69d8c10f1 100644 --- a/FIAT/quadrature_schemes.py +++ b/FIAT/quadrature_schemes.py @@ -62,14 +62,14 @@ def create_quadrature(ref_el, degree, scheme="default", entity=None): integration domain. If not provided, the domain is the entire cell. """ if entity is not None: - dim, entity_id = entity - sub_el = ref_el.construct_subelement(dim) + dimension, entity_id = entity + sub_el = ref_el.construct_subelement(dimension) Q_ref = create_quadrature(sub_el, degree, scheme=scheme) - return FacetQuadratureRule(ref_el, dim, entity_id, Q_ref) + return FacetQuadratureRule(ref_el, dimension, entity_id, Q_ref) if ref_el.is_macrocell(): - dim = ref_el.get_dimension() - sub_el = ref_el.construct_subelement(dim) + dimension = ref_el.get_dimension() + sub_el = ref_el.construct_subelement(dimension) Q_ref = create_quadrature(sub_el, degree, scheme=scheme) return MacroQuadratureRule(ref_el, Q_ref) From 053895f34581401fce8ce27daf2272eba6b6d237 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 27 Oct 2024 17:34:58 +0000 Subject: [PATCH 20/25] fix JM dof ordering --- FIAT/arnold_winther.py | 6 +++--- FIAT/hu_zhang.py | 3 ++- FIAT/johnson_mercier.py | 16 ++++++++-------- 3 files changed, 13 insertions(+), 12 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index aa2523e9a..fcb3b8708 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -9,7 +9,7 @@ # SPDX-License-Identifier: LGPL-3.0-or-later -from FIAT import finite_element, dual_set, expansions, polynomial_set +from FIAT import finite_element, dual_set, polynomial_set from FIAT.reference_element import TRIANGLE from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (ComponentPointEvaluation, @@ -109,9 +109,9 @@ def __init__(self, ref_el, degree=3): for phi in phis for i in range(sd) for j in range(i, sd)) # constraint dofs: moments of divergence against P_{k-1} \ P_{k-2} - dimPkm1 = expansions.polynomial_dimension(ref_el, degree-1) - dimPkm2 = expansions.polynomial_dimension(ref_el, degree-2) P = polynomial_set.ONPolynomialSet(ref_el, degree-1, shape=(sd,)) + dimPkm1 = P.expansion_set.get_num_members(degree-1) + dimPkm2 = P.expansion_set.get_num_members(degree-2) PH = P.take([i + j * dimPkm1 for j in range(sd) for i in range(dimPkm2, dimPkm1)]) phis = PH.tabulate(Q.get_points())[(0,)*sd] nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index ce7364e45..157f7ec3a 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -84,4 +84,5 @@ def __init__(self, ref_el, degree=3, variant=None): poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) dual = HuZhangDual(ref_el, degree, variant, qdegree) formdegree = ref_el.get_spatial_dimension() - 1 - super().__init__(poly_set, dual, degree, formdegree=formdegree, mapping="double contravariant piola") + mapping = "double contravariant piola" + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index de4683bf4..15473d796 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -38,13 +38,12 @@ def __init__(self, ref_complex, degree, variant=None): entity_ids[dim][facet].extend(range(cur, len(nodes))) cur = len(nodes) - if variant is None: - # Interior dofs: moments for each independent component - Q = create_quadrature(ref_complex, 2*degree-1) - P = polynomial_set.ONPolynomialSet(ref_el, degree-1) - phis = P.tabulate(Q.get_points())[(0,) * sd] - nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j)) - for j in range(sd) for i in range(j+1) for phi in phis) + # Interior dofs: moments for each independent component + Q = create_quadrature(ref_complex, 2*degree-1) + P = polynomial_set.ONPolynomialSet(ref_el, degree-1) + phis = P.tabulate(Q.get_points())[(0,) * sd] + nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j)) + for phi in phis for i in range(sd) for j in range(i, sd)) entity_ids[sd][0].extend(range(cur, len(nodes))) @@ -58,5 +57,6 @@ def __init__(self, ref_el, degree=1, variant=None): ref_complex = macro.AlfeldSplit(ref_el) poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) dual = JohnsonMercierDualSet(ref_complex, degree, variant=variant) + formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" - super().__init__(poly_set, dual, degree, mapping=mapping) + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) From ae5d811a05a67010491fd7a20529e57f6ab0bd3d Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Sun, 27 Oct 2024 23:34:30 +0000 Subject: [PATCH 21/25] new interior dofs --- FIAT/arnold_winther.py | 14 +++++++------- FIAT/expansions.py | 13 +++++++------ FIAT/functional.py | 2 +- FIAT/hu_zhang.py | 9 +++++---- FIAT/johnson_mercier.py | 24 +++++++++++------------- test/unit/test_awc.py | 28 +++++++++++++++++----------- test/unit/test_awnc.py | 22 ++++++++++++++-------- 7 files changed, 62 insertions(+), 50 deletions(-) diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index fcb3b8708..e7af81783 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -13,12 +13,11 @@ from FIAT.reference_element import TRIANGLE from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (ComponentPointEvaluation, - IntegralMoment, + TensorBidirectionalIntegralMoment, IntegralMomentOfTensorDivergence, IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment) - import numpy @@ -28,7 +27,6 @@ def __init__(self, ref_el, degree=2): raise ValueError("Nonconforming Arnold-Winther elements are only defined for degree 2.") top = ref_el.get_topology() sd = ref_el.get_spatial_dimension() - shp = (sd, sd) entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} nodes = [] @@ -45,9 +43,10 @@ def __init__(self, ref_el, degree=2): # internal dofs: constant moments of three unique components cur = len(nodes) + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) Q = create_quadrature(ref_el, degree) - phi = numpy.ones(Q.get_weights().shape) - nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + phi = numpy.full(Q.get_weights().shape, 1/ref_el.volume()) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) for i in range(sd) for j in range(i, sd)) entity_ids[2][0].extend(range(cur, len(nodes))) @@ -102,10 +101,11 @@ def __init__(self, ref_el, degree=3): entity_ids[1][entity].extend(range(cur, len(nodes))) # internal dofs: moments of unique components against P_{k-3} + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) Q = create_quadrature(ref_el, 2*(degree-1)) - P = polynomial_set.ONPolynomialSet(ref_el, degree-3) + P = polynomial_set.ONPolynomialSet(ref_el, degree-3, scale="L2 piola") phis = P.tabulate(Q.get_points())[(0,)*sd] - nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) for phi in phis for i in range(sd) for j in range(i, sd)) # constraint dofs: moments of divergence against P_{k-1} \ P_{k-2} diff --git a/FIAT/expansions.py b/FIAT/expansions.py index 7fa0ad887..28eb561d4 100644 --- a/FIAT/expansions.py +++ b/FIAT/expansions.py @@ -279,16 +279,19 @@ def __init__(self, ref_el, scale=None, variant=None): self._dmats_cache = {} self._cell_node_map_cache = {} - def get_scale(self, cell=0): + def get_scale(self, n, cell=0): scale = self.scale + sd = self.ref_el.get_spatial_dimension() if isinstance(scale, str): - sd = self.ref_el.get_spatial_dimension() vol = self.ref_el.volume_of_subcomplex(sd, cell) scale = scale.lower() if scale == "orthonormal": scale = math.sqrt(1.0 / vol) elif scale == "l2 piola": scale = 1.0 / vol + elif n == 0 and sd > 1 and len(self.affine_mappings) == 1: + # return 1 for n=0 to make regression tests pass + scale = 1 return scale def get_num_members(self, n): @@ -310,9 +313,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): ref_pts = numpy.add(numpy.dot(pts, A.T), b).T Jinv = A if direction is None else numpy.dot(A, direction)[:, None] sd = self.ref_el.get_spatial_dimension() - - # Always return 1 for n=0 to make regression tests pass - scale = 1.0 if n == 0 and len(self.affine_mappings) == 1 else self.get_scale(cell=cell) + scale = self.get_scale(n, cell=cell) phi = dubiner_recurrence(sd, n, lorder, ref_pts, Jinv, scale, variant=self.variant) if self.continuity == "C0": @@ -549,7 +550,7 @@ def _tabulate_on_cell(self, n, pts, order=0, cell=0, direction=None): Jinv = A[0, 0] if direction is None else numpy.dot(A, direction) xs = numpy.add(numpy.dot(pts, A.T), b) results = {} - scale = self.get_scale(cell=cell) * numpy.sqrt(2 * numpy.arange(n+1) + 1) + scale = self.get_scale(n, cell=cell) * 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/FIAT/functional.py b/FIAT/functional.py index 817701187..720412c5d 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -625,7 +625,7 @@ def __init__(self, ref_el, v, w, pt): super().__init__(ref_el, shp, pt_dict, {}, "PointwiseInnerProductEval") -class TensorBidirectionalMomentInnerProductEvaluation(FrobeniusIntegralMoment): +class TensorBidirectionalIntegralMoment(FrobeniusIntegralMoment): r""" This is a functional on symmetric 2-tensor fields. Let u be such a field, f a function tabulated at points, and v,w be vectors. This implements the evaluation diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py index 157f7ec3a..5068a6b39 100644 --- a/FIAT/hu_zhang.py +++ b/FIAT/hu_zhang.py @@ -14,7 +14,7 @@ from FIAT.quadrature_schemes import create_quadrature from FIAT.functional import (ComponentPointEvaluation, PointwiseInnerProductEvaluation, - IntegralMoment, + TensorBidirectionalIntegralMoment, IntegralLegendreNormalNormalMoment, IntegralLegendreNormalTangentialMoment) @@ -63,10 +63,11 @@ def __init__(self, ref_el, degree, variant, qdegree): elif variant == "integral": # Moments of unique components against a basis for P_{k-2} - Q = create_quadrature(ref_el, qdegree + degree-2) - P = polynomial_set.ONPolynomialSet(ref_el, degree-2) + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) + Q = create_quadrature(ref_el, 2*degree-2) + P = polynomial_set.ONPolynomialSet(ref_el, degree-2, scale="L2 piola") phis = P.tabulate(Q.get_points())[(0,)*sd] - nodes.extend(IntegralMoment(ref_el, Q, phi, (i, j), shp) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) for phi in phis for i in range(sd) for j in range(i, sd)) entity_ids[2][0].extend(range(cur, len(nodes))) diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 15473d796..a36907212 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -1,5 +1,5 @@ from FIAT import finite_element, dual_set, macro, polynomial_set -from FIAT.functional import IntegralMoment, FrobeniusIntegralMoment +from FIAT.functional import TensorBidirectionalIntegralMoment from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature import numpy @@ -18,31 +18,29 @@ def __init__(self, ref_complex, degree, variant=None): nodes = [] # Face dofs: bidirectional (nn and nt) Legendre moments - R = numpy.array([[0, 1], [-1, 0]]) dim = sd - 1 + R = numpy.array([[0, 1], [-1, 0]]) ref_facet = ref_el.construct_subelement(dim) Qref = create_quadrature(ref_facet, 2*degree) P = polynomial_set.ONPolynomialSet(ref_facet, degree) phis = P.tabulate(Qref.get_points())[(0,) * dim] - for facet in sorted(top[dim]): + for f in sorted(top[dim]): cur = len(nodes) - Q = FacetQuadratureRule(ref_el, dim, facet, Qref) - thats = ref_el.compute_tangents(dim, facet) + Q = FacetQuadratureRule(ref_el, dim, f, Qref) + thats = ref_el.compute_tangents(dim, f) nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) normal = nhat / Q.jacobian_determinant() - - uvecs = (nhat, *thats) - comps = [numpy.outer(normal, uvec) for uvec in uvecs] - nodes.extend(FrobeniusIntegralMoment(ref_el, Q, comp[:, :, None] * phi[None, None, :]) - for phi in phis for comp in comps) - entity_ids[dim][facet].extend(range(cur, len(nodes))) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, normal, comp, Q, phi) + for phi in phis for comp in (nhat, *thats)) + entity_ids[dim][f].extend(range(cur, len(nodes))) cur = len(nodes) # Interior dofs: moments for each independent component + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) Q = create_quadrature(ref_complex, 2*degree-1) - P = polynomial_set.ONPolynomialSet(ref_el, degree-1) + P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola") phis = P.tabulate(Q.get_points())[(0,) * sd] - nodes.extend(IntegralMoment(ref_el, Q, phi, comp=(i, j)) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n[i+1], n[j+1], Q, phi) for phi in phis for i in range(sd) for j in range(i, sd)) entity_ids[sd][0].extend(range(cur, len(nodes))) diff --git a/test/unit/test_awc.py b/test/unit/test_awc.py index c5387803a..b96bc86be 100644 --- a/test/unit/test_awc.py +++ b/test/unit/test_awc.py @@ -1,5 +1,5 @@ import numpy as np -from FIAT import ufc_simplex, ArnoldWinther, make_quadrature, expansions +from FIAT import ufc_simplex, ArnoldWinther, create_quadrature, expansions def test_dofs(): @@ -20,7 +20,7 @@ def test_dofs(): assert np.allclose(vert_vals[3*i+j, :, :, (i+k) % 3], np.zeros((2, 2))) # check edge moments - Qline = make_quadrature(line, 6) + Qline = create_quadrature(line, 6) linebfs = expansions.LineExpansionSet(line) linevals = linebfs.tabulate(1, Qline.pts) @@ -73,15 +73,21 @@ def test_dofs(): assert np.allclose(ntmoments[bf, :], np.zeros(2)) # check internal dofs - Q = make_quadrature(T, 6) + ns = list(map(T.compute_scaled_normal, range(3))) + Q = create_quadrature(T, 3) qpvals = AW.tabulate(0, Q.pts)[(0, 0)] - const_moms = qpvals @ Q.wts - assert np.allclose(const_moms[:21], np.zeros((21, 2, 2))) - assert np.allclose(const_moms[24:], np.zeros((6, 2, 2))) - assert np.allclose(const_moms[21:24, 0, 0], np.asarray([1, 0, 0])) - assert np.allclose(const_moms[21:24, 0, 1], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[21:24, 1, 0], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[21:24, 1, 1], np.asarray([0, 0, 1])) + const_moms = qpvals @ Q.wts / T.volume() + nn_moms = const_moms.copy() + for j in range(2): + for i in range(2): + comp = np.outer(ns[i+1], ns[j+1]) + nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1))) + assert np.allclose(nn_moms[:21], np.zeros((21, 2, 2))) + assert np.allclose(nn_moms[24:], np.zeros((6, 2, 2))) + assert np.allclose(nn_moms[21:24, 0, 0], np.asarray([1, 0, 0])) + assert np.allclose(nn_moms[21:24, 0, 1], np.asarray([0, 1, 0])) + assert np.allclose(nn_moms[21:24, 1, 0], np.asarray([0, 1, 0])) + assert np.allclose(nn_moms[21:24, 1, 1], np.asarray([0, 0, 1])) def frob(a, b): @@ -94,7 +100,7 @@ def test_projection(): AW = ArnoldWinther(T, 3) - Q = make_quadrature(T, 4) + Q = create_quadrature(T, 6) qpts = np.asarray(Q.pts) qwts = np.asarray(Q.wts) nqp = len(Q.wts) diff --git a/test/unit/test_awnc.py b/test/unit/test_awnc.py index f7a28cac8..c770791ee 100644 --- a/test/unit/test_awnc.py +++ b/test/unit/test_awnc.py @@ -5,7 +5,7 @@ def test_dofs(): line = ufc_simplex(1) T = ufc_simplex(2) - T.vertices = np.asarray([(0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)]) + T.vertices = ((0.0, 0.0), (1.0, 0.25), (-0.75, 1.1)) AW = ArnoldWintherNC(T, 2) Qline = make_quadrature(line, 6) @@ -61,12 +61,18 @@ def test_dofs(): assert np.allclose(ntmoments[bf, :], np.zeros(2), atol=1.e-7) # check internal dofs + ns = list(map(T.compute_scaled_normal, range(3))) Q = make_quadrature(T, 6) qpvals = AW.tabulate(0, Q.pts)[(0, 0)] - const_moms = qpvals @ Q.wts - assert np.allclose(const_moms[:12], np.zeros((12, 2, 2))) - assert np.allclose(const_moms[15:], np.zeros((3, 2, 2))) - assert np.allclose(const_moms[12:15, 0, 0], np.asarray([1, 0, 0])) - assert np.allclose(const_moms[12:15, 0, 1], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[12:15, 1, 0], np.asarray([0, 1, 0])) - assert np.allclose(const_moms[12:15, 1, 1], np.asarray([0, 0, 1])) + const_moms = qpvals @ Q.wts / T.volume() + nn_moms = const_moms.copy() + for j in range(2): + for i in range(2): + comp = np.outer(ns[i+1], ns[j+1]) + nn_moms[:, i, j] = np.tensordot(const_moms, comp, ((1, 2), (0, 1))) + assert np.allclose(nn_moms[:12], np.zeros((12, 2, 2))) + assert np.allclose(nn_moms[15:], np.zeros((3, 2, 2))) + assert np.allclose(nn_moms[12:15, 0, 0], np.asarray([1, 0, 0])) + assert np.allclose(nn_moms[12:15, 0, 1], np.asarray([0, 1, 0])) + assert np.allclose(nn_moms[12:15, 1, 0], np.asarray([0, 1, 0])) + assert np.allclose(nn_moms[12:15, 1, 1], np.asarray([0, 0, 1])) From a46bd3640463394cd5acc702ccb8e102a3808bf4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 29 Oct 2024 17:47:20 +0000 Subject: [PATCH 22/25] address review comments --- FIAT/mardal_tai_winther.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 354a71601..bb9b652d0 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -66,9 +66,9 @@ def __init__(self, cell, degree): dofs.extend(_dofs) for entity_id in range(3): - dof_ids[1][entity_id] = dof_ids[1][entity_id] + _edge_dof_ids[entity_id] + dof_ids[1][entity_id].extend(_edge_dof_ids[entity_id]) - dof_ids[2][0] = dof_ids[2][0] + _cell_dof_ids + dof_ids[2][0].extend(_cell_dof_ids) super(MardalTaiWintherDual, self).__init__(dofs, cell, dof_ids) @@ -118,18 +118,21 @@ def _generate_constraint_dofs(cell, degree, offset): as described in the FIAT paper. """ dofs = [] - edge_dof_ids = {} + + start_order = 2 + stop_order = 3 + qdegree = degree + stop_order for entity_id in range(3): - dofs.append(IntegralLegendreNormalMoment(cell, entity_id, 2, degree+3)) - dofs.append(IntegralLegendreNormalMoment(cell, entity_id, 3, degree+3)) + cur = len(dofs) + dofs.extend(IntegralLegendreNormalMoment(cell, entity_id, order, qdegree) + for order in range(start_order, stop_order+1)) - edge_dof_ids[entity_id] = [offset, offset+1] - offset += 2 + edge_dof_ids[entity_id] = list(range(offset+cur, offset+len(dofs))) - cell_dofs = DivergenceDubinerMoments(cell, 1, 2, 6) - dofs.extend(cell_dofs) - cell_dof_ids = list(range(offset, offset+len(cell_dofs))) + cur = len(dofs) + dofs.extend(DivergenceDubinerMoments(cell, start_order-1, stop_order-1, degree)) + cell_dof_ids = list(range(offset+cur, offset+len(dofs))) return (dofs, edge_dof_ids, cell_dof_ids) From 98b86b32184eb413f631d30cb80c31c5587c1063 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 1 Nov 2024 14:30:05 +0000 Subject: [PATCH 23/25] cleanup MTW --- FIAT/mardal_tai_winther.py | 176 +++++++++++++------------------------ 1 file changed, 62 insertions(+), 114 deletions(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index bb9b652d0..9188705b6 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -17,140 +17,88 @@ from FIAT.quadrature_schemes import create_quadrature -def DivergenceDubinerMoments(cell, start_deg, stop_deg, comp_deg): - sd = cell.get_spatial_dimension() - P = polynomial_set.ONPolynomialSet(cell, stop_deg) - Q = create_quadrature(cell, comp_deg + stop_deg) +def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg): + sd = ref_el.get_spatial_dimension() + P = polynomial_set.ONPolynomialSet(ref_el, stop_deg) + Q = create_quadrature(ref_el, comp_deg + stop_deg) - dim0 = expansions.polynomial_dimension(cell, start_deg-1) - dim1 = expansions.polynomial_dimension(cell, stop_deg) + dim0 = expansions.polynomial_dimension(ref_el, start_deg-1) + dim1 = expansions.polynomial_dimension(ref_el, stop_deg) indices = list(range(dim0, dim1)) phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd] - return [IntegralMomentOfDivergence(cell, Q, phi) for phi in phis] + return [IntegralMomentOfDivergence(ref_el, Q, phi) for phi in phis] class MardalTaiWintherDual(dual_set.DualSet): """Degrees of freedom for Mardal-Tai-Winther elements.""" - def __init__(self, cell, degree): - dim = cell.get_spatial_dimension() - if not dim == 2: - raise ValueError("Mardal-Tai-Winther elements are only" - "defined in dimension 2.") - - if not degree == 3: - raise ValueError("Mardal-Tai-Winther elements are only defined" - "for degree 3.") - - # construct the degrees of freedoms - dofs = [] # list of functionals - - # dof_ids[i][j] contains the indices of dofs that are associated with - # entity j in dim i - dof_ids = {} - - # no vertex dof - dof_ids[0] = {i: [] for i in range(dim + 1)} - - # edge dofs - (_dofs, _dof_ids) = self._generate_edge_dofs(cell, degree) - dofs.extend(_dofs) - dof_ids[1] = _dof_ids - - # no cell dofs - dof_ids[2] = {} - dof_ids[2][0] = [] - - # extra dofs for enforcing div(v) constant over the cell and - # v.n linear on edges - (_dofs, _edge_dof_ids, _cell_dof_ids) = self._generate_constraint_dofs(cell, degree, len(dofs)) - dofs.extend(_dofs) - - for entity_id in range(3): - dof_ids[1][entity_id].extend(_edge_dof_ids[entity_id]) - - dof_ids[2][0].extend(_cell_dof_ids) - - super(MardalTaiWintherDual, self).__init__(dofs, cell, dof_ids) - - @staticmethod - def _generate_edge_dofs(cell, degree): - """Generate dofs on edges. - On each edge, let n be its normal. We need to integrate - u.n and u.t against the first Legendre polynomial (constant) - and u.n against the second (linear). - """ - dofs = [] - dof_ids = {} - offset = 0 - sd = 2 - - facet = cell.get_facet_element() - # Facet nodes are \int_F v\cdot n p ds where p \in P_{q-1} + def __init__(self, ref_el, degree): + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + + if sd != 2: + raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2.") + + if degree != 3: + raise ValueError("Mardal-Tai-Winther elements are only defined for degree 3.") + + entity_ids = {dim: {entity: [] for entity in top[dim]} for dim in top} + nodes = [] + + # no vertex dofs + + # On each facet, let n be its normal. We need to integrate + # u.n and u.t against the first Legendre polynomial (constant) + # and u.n against the second (linear). + facet = ref_el.get_facet_element() + # Facet nodes are \int_F v.n p ds where p \in P_{q-1} # degree is q - 1 Q = create_quadrature(facet, degree+1) Pq = polynomial_set.ONPolynomialSet(facet, 1) phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)] - for f in range(3): - dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phis[0], f)) - dofs.append(IntegralMomentOfTangentialEvaluation(cell, Q, phis[0], f)) - dofs.append(IntegralMomentOfNormalEvaluation(cell, Q, phis[1], f)) - - num_new_dofs = 3 - dof_ids[f] = list(range(offset, offset + num_new_dofs)) - offset += num_new_dofs - - return (dofs, dof_ids) - - @staticmethod - def _generate_constraint_dofs(cell, degree, offset): - """ - Generate constraint dofs on the cell and edges - * div(v) must be constant on the cell. Since v is a cubic and - div(v) is quadratic, we need the integral of div(v) against the - linear and quadratic Dubiner polynomials to vanish. - There are two linear and three quadratics, so these are five - constraints - * v.n must be linear on each edge. Since v.n is cubic, we need - the integral of v.n against the cubic and quadratic Legendre - polynomial to vanish on each edge. - - So we introduce functionals whose kernel describes this property, - as described in the FIAT paper. - """ - dofs = [] - edge_dof_ids = {} - + for f in sorted(top[sd-1]): + cur = len(nodes) + nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[0], f)) + nodes.append(IntegralMomentOfTangentialEvaluation(ref_el, Q, phis[0], f)) + nodes.append(IntegralMomentOfNormalEvaluation(ref_el, Q, phis[1], f)) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + # Generate constraint nodes on the cell and facets + # * div(v) must be constant on the cell. Since v is a cubic and + # div(v) is quadratic, we need the integral of div(v) against the + # linear and quadratic Dubiner polynomials to vanish. + # There are two linear and three quadratics, so these are five + # constraints + # * v.n must be linear on each facet. Since v.n is cubic, we need + # the integral of v.n against the cubic and quadratic Legendre + # polynomial to vanish on each facet. + + # So we introduce functionals whose kernel describes this property, + # as described in the FIAT paper. start_order = 2 stop_order = 3 qdegree = degree + stop_order - for entity_id in range(3): - cur = len(dofs) - dofs.extend(IntegralLegendreNormalMoment(cell, entity_id, order, qdegree) - for order in range(start_order, stop_order+1)) - - edge_dof_ids[entity_id] = list(range(offset+cur, offset+len(dofs))) + for f in sorted(top[sd-1]): + cur = len(nodes) + nodes.extend(IntegralLegendreNormalMoment(ref_el, f, order, qdegree) + for order in range(start_order, stop_order+1)) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) - cur = len(dofs) - dofs.extend(DivergenceDubinerMoments(cell, start_order-1, stop_order-1, degree)) - cell_dof_ids = list(range(offset+cur, offset+len(dofs))) + cur = len(nodes) + nodes.extend(DivergenceDubinerMoments(ref_el, start_order-1, stop_order-1, degree)) + entity_ids[sd][0].extend(range(cur, len(nodes))) - return (dofs, edge_dof_ids, cell_dof_ids) + super().__init__(nodes, ref_el, entity_ids) class MardalTaiWinther(finite_element.CiarletElement): """The definition of the Mardal-Tai-Winther element. """ - def __init__(self, cell, degree=3): + def __init__(self, ref_el, degree=3): + sd = ref_el.get_spatial_dimension() assert degree == 3, "Only defined for degree 3" - assert cell.get_spatial_dimension() == 2, "Only defined for dimension 2" - # polynomial space - Ps = polynomial_set.ONPolynomialSet(cell, degree, (2,)) - - # degrees of freedom - Ls = MardalTaiWintherDual(cell, degree) - - # mapping under affine transformation + assert sd == 2, "Only defined for dimension 2" + poly_set = polynomial_set.ONPolynomialSet(ref_el, degree, (sd,)) + dual = MardalTaiWintherDual(ref_el, degree) + formdegree = sd-1 mapping = "contravariant piola" - - super(MardalTaiWinther, self).__init__(Ps, Ls, degree, - mapping=mapping) + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) From d8c61e18cfe7471ef65cbdddb2505d91265fd091 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Fri, 1 Nov 2024 15:04:31 +0000 Subject: [PATCH 24/25] This is a FunctionalFactory --- FIAT/mardal_tai_winther.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 9188705b6..2c623c64f 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -26,7 +26,8 @@ def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg): dim1 = expansions.polynomial_dimension(ref_el, stop_deg) indices = list(range(dim0, dim1)) phis = P.take(indices).tabulate(Q.get_points())[(0,)*sd] - return [IntegralMomentOfDivergence(ref_el, Q, phi) for phi in phis] + for phi in phis: + yield IntegralMomentOfDivergence(ref_el, Q, phi) class MardalTaiWintherDual(dual_set.DualSet): From d25817b50d849cc7e48481932f179b349b717abd Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Wed, 6 Nov 2024 23:25:17 +0000 Subject: [PATCH 25/25] Implement a traceless H(curl div) element (#69) * TracelessTensorPolynomialSet * Add GLS H(curl div) element --- FIAT/__init__.py | 4 + FIAT/functional.py | 6 +- FIAT/gopalakrishnan_lederer_schoberl.py | 88 +++++++++++++++++++ FIAT/guzman_neilan.py | 5 +- FIAT/polynomial_set.py | 67 +++++++++----- test/unit/test_fiat.py | 14 +++ .../test_gopalakrishnan_lederer_schoberl.py | 77 ++++++++++++++++ 7 files changed, 234 insertions(+), 27 deletions(-) create mode 100644 FIAT/gopalakrishnan_lederer_schoberl.py create mode 100644 test/unit/test_gopalakrishnan_lederer_schoberl.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 8c9fcc1e5..26c467a99 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -48,6 +48,8 @@ from FIAT.raviart_thomas import RaviartThomas from FIAT.crouzeix_raviart import CrouzeixRaviart from FIAT.regge import Regge +from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind +from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson from FIAT.arnold_winther import ArnoldWinther from FIAT.arnold_winther import ArnoldWintherNC @@ -108,6 +110,8 @@ "BrokenElement": DiscontinuousElement, "HDiv Trace": HDivTrace, "Hellan-Herrmann-Johnson": HellanHerrmannJohnson, + "Gopalakrishnan-Lederer-Schoberl 1st kind": GopalakrishnanLedererSchoberlFirstKind, + "Gopalakrishnan-Lederer-Schoberl 2nd kind": GopalakrishnanLedererSchoberlSecondKind, "Conforming Arnold-Winther": ArnoldWinther, "Nonconforming Arnold-Winther": ArnoldWintherNC, "Hu-Zhang": HuZhang, diff --git a/FIAT/functional.py b/FIAT/functional.py index 720412c5d..6a831a87c 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -631,13 +631,13 @@ class TensorBidirectionalIntegralMoment(FrobeniusIntegralMoment): field, f a function tabulated at points, and v,w be vectors. This implements the evaluation \int v^T u(x) w f(x). Clearly v^iu_{ij}w^j = u_{ij}v^iw^j. Thus the value can be computed - from the Frobenius inner product of u with wv^T. This gives the + from the Frobenius inner product of u with vw^T. This gives the correct weights. """ def __init__(self, ref_el, v, w, Q, f_at_qpts): - wvT = numpy.outer(w, v) - F_at_qpts = numpy.multiply(wvT[..., None], f_at_qpts) + vwT = numpy.outer(v, w) + F_at_qpts = numpy.multiply(vwT[..., None], f_at_qpts) super().__init__(ref_el, Q, F_at_qpts, "TensorBidirectionalMomentInnerProductEvaluation") diff --git a/FIAT/gopalakrishnan_lederer_schoberl.py b/FIAT/gopalakrishnan_lederer_schoberl.py new file mode 100644 index 000000000..656cee6bb --- /dev/null +++ b/FIAT/gopalakrishnan_lederer_schoberl.py @@ -0,0 +1,88 @@ +from FIAT import finite_element, dual_set, polynomial_set, expansions +from FIAT.functional import TensorBidirectionalIntegralMoment as BidirectionalMoment +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule +from FIAT.restricted import RestrictedElement + + +class GLSDual(dual_set.DualSet): + def __init__(self, ref_el, degree): + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + nodes = [] + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + + # Face dofs: moments of normal-tangential components against a basis for Pk + ref_facet = ref_el.construct_subelement(sd-1) + Qref = create_quadrature(ref_facet, 2*degree) + P = polynomial_set.ONPolynomialSet(ref_facet, degree) + phis = P.tabulate(Qref.get_points())[(0,) * (sd-1)] + for f in sorted(top[sd-1]): + cur = len(nodes) + Q = FacetQuadratureRule(ref_el, sd-1, f, Qref) + Jdet = Q.jacobian_determinant() + normal = ref_el.compute_scaled_normal(f) + tangents = ref_el.compute_tangents(sd-1, f) + n = normal / Jdet + nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi) + for phi in phis for t in tangents) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + # Interior dofs: moments of normal-tangential components against a basis for P_{k-1} + if degree > 0: + cur = len(nodes) + Q = create_quadrature(ref_el, 2*degree-1) + P = polynomial_set.ONPolynomialSet(ref_el, degree-1, scale="L2 piola") + phis = P.tabulate(Q.get_points())[(0,) * sd] + for f in sorted(top[sd-1]): + n = ref_el.compute_scaled_normal(f) + tangents = ref_el.compute_tangents(sd-1, f) + nodes.extend(BidirectionalMoment(ref_el, t, n, Q, phi) + for phi in phis for t in tangents) + entity_ids[sd][0].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) + + +class GopalakrishnanLedererSchoberlSecondKind(finite_element.CiarletElement): + """The GLS element used for the Mass-Conserving mixed Stress (MCS) + formulation for Stokes flow with weakly imposed stress symmetry. + + GLS^2(k) is the space of trace-free polynomials of degree k with + continuous normal-tangential components. + + Reference: https://doi.org/10.1137/19M1248960 + + Notes + ----- + This element does not include the bubbles required for inf-sup stability of + the weak symmetry constraint. + + """ + def __init__(self, ref_el, degree): + poly_set = polynomial_set.TracelessTensorPolynomialSet(ref_el, degree) + dual = GLSDual(ref_el, degree) + sd = ref_el.get_spatial_dimension() + formdegree = (1, sd-1) + mapping = "covariant contravariant piola" + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) + + +def GopalakrishnanLedererSchoberlFirstKind(ref_el, degree): + """The GLS element used for the Mass-Conserving mixed Stress (MCS) + formulation for Stokes flow. + + GLS^1(k) is the space of trace-free polynomials of degree k with + continuous normal-tangential components of degree k-1. + + Reference: https://doi.org/10.1093/imanum/drz022 + """ + fe = GopalakrishnanLedererSchoberlSecondKind(ref_el, degree) + entity_dofs = fe.entity_dofs() + sd = ref_el.get_spatial_dimension() + dimPkm1 = (sd-1)*expansions.polynomial_dimension(ref_el.construct_subelement(sd-1), degree-1) + indices = [] + for f in entity_dofs[sd-1]: + indices.extend(entity_dofs[sd-1][f][:dimPkm1]) + indices.extend(entity_dofs[sd][0]) + return RestrictedElement(fe, indices=indices) diff --git a/FIAT/guzman_neilan.py b/FIAT/guzman_neilan.py index ff853e583..970dbdf89 100644 --- a/FIAT/guzman_neilan.py +++ b/FIAT/guzman_neilan.py @@ -109,15 +109,12 @@ def GuzmanNeilanH1div(ref_el, degree=2, reduced=False): """ order = 0 AS = AlfeldSorokina(ref_el, 2) - if reduced: + if reduced or ref_el.get_spatial_dimension() <= 2: order = 1 # Only extract the div bubbles div_nodes = [i for i, node in enumerate(AS.dual_basis()) if len(node.deriv_dict) > 0] AS = RestrictedElement(AS, indices=div_nodes) - elif ref_el.get_spatial_dimension() <= 2: - # Quadratic bubbles are already included in 2D - return AS GN = GuzmanNeilanH1(ref_el, order=order) return NodalEnrichedElement(AS, GN) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 579ec3614..6210ef0b4 100644 --- a/FIAT/polynomial_set.py +++ b/FIAT/polynomial_set.py @@ -125,16 +125,13 @@ def __init__(self, ref_el, degree, shape=tuple(), **kwargs): if shape == tuple(): coeffs = numpy.eye(num_members) else: - coeffs_shape = (num_members, *shape, num_exp_functions) - coeffs = numpy.zeros(coeffs_shape, "d") - # use functional's index_iterator function - cur_bf = 0 + coeffs = numpy.zeros((num_members, *shape, num_exp_functions)) + cur = 0 + exp_bf = range(num_exp_functions) for idx in index_iterator(shape): - n = expansion_set.get_num_members(embedded_degree) - for exp_bf in range(n): - cur_idx = (cur_bf, *idx, exp_bf) - coeffs[cur_idx] = 1.0 - cur_bf += 1 + cur_bf = range(cur, cur+num_exp_functions) + coeffs[(cur_bf, *idx, exp_bf)] = 1.0 + cur += num_exp_functions super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) @@ -209,19 +206,49 @@ def __init__(self, ref_el, degree, size=None, **kwargs): embedded_degree = degree # set up coefficients for symmetric tensors - coeffs_shape = (num_members, *shape, num_exp_functions) - coeffs = numpy.zeros(coeffs_shape, "d") - cur_bf = 0 + coeffs = numpy.zeros((num_members, *shape, num_exp_functions)) + cur = 0 + exp_bf = range(num_exp_functions) for i, j in index_iterator(shape): + if i > j: + continue + cur_bf = range(cur, cur+num_exp_functions) + coeffs[cur_bf, i, j, exp_bf] = 1.0 + coeffs[cur_bf, j, i, exp_bf] = 1.0 + cur += num_exp_functions + + super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) + + +class TracelessTensorPolynomialSet(PolynomialSet): + """Constructs an orthonormal basis for traceless-tensor-valued + polynomials on a reference element. + """ + def __init__(self, ref_el, degree, size=None, **kwargs): + expansion_set = expansions.ExpansionSet(ref_el, **kwargs) + + sd = ref_el.get_spatial_dimension() + if size is None: + size = sd + + shape = (size, size) + num_exp_functions = expansion_set.get_num_members(degree) + num_components = size * size - 1 + num_members = num_components * num_exp_functions + embedded_degree = degree + + # set up coefficients for traceless tensors + coeffs = numpy.zeros((num_members, *shape, num_exp_functions)) + cur = 0 + exp_bf = range(num_exp_functions) + for i, j in index_iterator(shape): + if i == size-1 and j == size-1: + continue + cur_bf = range(cur, cur+num_exp_functions) + coeffs[cur_bf, i, j, exp_bf] = 1.0 if i == j: - for exp_bf in range(num_exp_functions): - coeffs[cur_bf, i, j, exp_bf] = 1.0 - cur_bf += 1 - elif i < j: - for exp_bf in range(num_exp_functions): - coeffs[cur_bf, i, j, exp_bf] = 1.0 - coeffs[cur_bf, j, i, exp_bf] = 1.0 - cur_bf += 1 + coeffs[cur_bf, -1, -1, exp_bf] = -1.0 + cur += num_exp_functions super().__init__(ref_el, degree, embedded_degree, expansion_set, coeffs) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index 158070e35..806e74911 100644 --- a/test/unit/test_fiat.py +++ b/test/unit/test_fiat.py @@ -35,6 +35,8 @@ from FIAT.regge import Regge # noqa: F401 from FIAT.hdiv_trace import HDivTrace, map_to_reference_facet # noqa: F401 from FIAT.hellan_herrmann_johnson import HellanHerrmannJohnson # noqa: F401 +from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlFirstKind # noqa: F401 +from FIAT.gopalakrishnan_lederer_schoberl import GopalakrishnanLedererSchoberlSecondKind # noqa: F401 from FIAT.brezzi_douglas_fortin_marini import BrezziDouglasFortinMarini # noqa: F401 from FIAT.gauss_legendre import GaussLegendre # noqa: F401 from FIAT.gauss_lobatto_legendre import GaussLobattoLegendre # noqa: F401 @@ -273,6 +275,18 @@ def __init__(self, a, b): "HellanHerrmannJohnson(S, 2)", "HellanHerrmannJohnson(T, 1, variant='point')", "HellanHerrmannJohnson(S, 1, variant='point')", + "GopalakrishnanLedererSchoberlFirstKind(T, 1)", + "GopalakrishnanLedererSchoberlFirstKind(T, 2)", + "GopalakrishnanLedererSchoberlFirstKind(T, 3)", + "GopalakrishnanLedererSchoberlFirstKind(S, 1)", + "GopalakrishnanLedererSchoberlFirstKind(S, 2)", + "GopalakrishnanLedererSchoberlFirstKind(S, 3)", + "GopalakrishnanLedererSchoberlSecondKind(T, 0)", + "GopalakrishnanLedererSchoberlSecondKind(T, 1)", + "GopalakrishnanLedererSchoberlSecondKind(T, 2)", + "GopalakrishnanLedererSchoberlSecondKind(S, 0)", + "GopalakrishnanLedererSchoberlSecondKind(S, 1)", + "GopalakrishnanLedererSchoberlSecondKind(S, 2)", "BrezziDouglasFortinMarini(T, 2)", "GaussLegendre(I, 0)", "GaussLegendre(I, 1)", diff --git a/test/unit/test_gopalakrishnan_lederer_schoberl.py b/test/unit/test_gopalakrishnan_lederer_schoberl.py new file mode 100644 index 000000000..f78c8da34 --- /dev/null +++ b/test/unit/test_gopalakrishnan_lederer_schoberl.py @@ -0,0 +1,77 @@ +import pytest +import numpy + +from FIAT import (GopalakrishnanLedererSchoberlFirstKind, + GopalakrishnanLedererSchoberlSecondKind) +from FIAT.reference_element import ufc_simplex +from FIAT.expansions import polynomial_dimension +from FIAT.polynomial_set import ONPolynomialSet +from FIAT.quadrature_schemes import create_quadrature +from FIAT.quadrature import FacetQuadratureRule + + +@pytest.fixture(params=("T", "S")) +def cell(request): + dim = {"I": 1, "T": 2, "S": 3}[request.param] + return ufc_simplex(dim) + + +@pytest.mark.parametrize("degree", (1, 2, 3)) +@pytest.mark.parametrize("kind", (1, 2)) +def test_gls_bubbles(kind, cell, degree): + if kind == 1: + element = GopalakrishnanLedererSchoberlFirstKind + else: + element = GopalakrishnanLedererSchoberlSecondKind + fe = element(cell, degree) + sd = cell.get_spatial_dimension() + facet_el = cell.construct_subelement(sd-1) + poly_set = fe.get_nodal_basis() + + # test dimension of constrained space + dimPkm1 = polynomial_dimension(facet_el, degree-1) + dimPkp1 = polynomial_dimension(facet_el, degree+1) + dimPk = polynomial_dimension(facet_el, degree) + if kind == 1: + constraints = dimPk - dimPkm1 + else: + constraints = 0 + expected = (sd**2-1)*(polynomial_dimension(cell, degree) - constraints) + assert poly_set.get_num_members() == expected + + # test dimension of the bubbles + entity_dofs = fe.entity_dofs() + bubbles = poly_set.take(entity_dofs[sd][0]) + expected = (sd**2-1)*polynomial_dimension(cell, degree-1) + assert bubbles.get_num_members() == expected + + top = cell.get_topology() + Qref = create_quadrature(facet_el, 2*degree+1) + Pk = ONPolynomialSet(facet_el, degree+1) + if kind == 1: + start, stop = dimPkm1, dimPkp1 + else: + start, stop = dimPk, dimPkp1 + PkH = Pk.take(list(range(start, stop))) + PkH_at_qpts = PkH.tabulate(Qref.get_points())[(0,)*(sd-1)] + weights = numpy.transpose(numpy.multiply(PkH_at_qpts, Qref.get_weights())) + for facet in top[sd-1]: + n = cell.compute_scaled_normal(facet) + rts = cell.compute_tangents(sd-1, facet) + Q = FacetQuadratureRule(cell, sd-1, facet, Qref) + qpts, qwts = Q.get_points(), Q.get_weights() + + # test the degree of normal-tangential components + phi_at_pts = fe.tabulate(0, qpts)[(0,) * sd] + for t in rts: + nt = numpy.outer(t, n) + phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2))) + assert numpy.allclose(numpy.dot(phi_nt, weights), 0) + + # test the support of the normal-tangential bubble + phi_at_pts = bubbles.tabulate(qpts)[(0,) * sd] + for t in rts: + nt = numpy.outer(t, n) + phi_nt = numpy.tensordot(nt, phi_at_pts, axes=((0, 1), (1, 2))) + norms = numpy.dot(phi_nt**2, qwts) + assert numpy.allclose(norms, 0)