From b06385456598b185e253e6e1d6eb7719407e29ca Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Tue, 19 Nov 2024 17:09:52 +0000 Subject: [PATCH] Integral variant for Regge/HHJ (#96) * docstring * cleanup test_fiat.py * Regge integral dofs * HHJ integral dofs in 2d/3d * support variant="integral(q)" * BidirectionalMoment dofs * address review comments * Cleanup H(div, S) elements (#98) * implementation of Hu-Zhang element --------- Co-authored-by: Francis Aznaran --- FIAT/__init__.py | 33 +-- FIAT/arnold_winther.py | 228 ++++++++---------- FIAT/expansions.py | 13 +- FIAT/functional.py | 192 +++++---------- FIAT/gopalakrishnan_lederer_schoberl.py | 88 +++++++ FIAT/guzman_neilan.py | 5 +- FIAT/hellan_herrmann_johnson.py | 131 +++++----- FIAT/hu_zhang.py | 89 +++++++ FIAT/johnson_mercier.py | 36 ++- FIAT/mardal_tai_winther.py | 210 ++++++---------- FIAT/polynomial_set.py | 67 +++-- FIAT/quadrature_schemes.py | 40 ++- FIAT/reference_element.py | 2 +- FIAT/regge.py | 78 +++--- test/regression/test_regression.py | 5 +- test/unit/test_awc.py | 32 ++- test/unit/test_awnc.py | 22 +- test/unit/test_fiat.py | 123 +++------- .../test_gopalakrishnan_lederer_schoberl.py | 77 ++++++ test/unit/test_polynomial.py | 109 +++++++++ 20 files changed, 909 insertions(+), 671 deletions(-) create mode 100644 FIAT/gopalakrishnan_lederer_schoberl.py create mode 100644 FIAT/hu_zhang.py create mode 100644 test/unit/test_gopalakrishnan_lederer_schoberl.py create mode 100644 test/unit/test_polynomial.py diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 117606fd7..26c467a99 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -2,10 +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.finite_element import FiniteElement, CiarletElement # 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 # 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 @@ -17,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 @@ -42,9 +48,12 @@ 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 +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 @@ -52,20 +61,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 -# 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 supported_elements = {"Argyris": Argyris, "Bell": Bell, @@ -112,8 +110,11 @@ "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, "Mardal-Tai-Winther": MardalTaiWinther} # List of extra elements diff --git a/FIAT/arnold_winther.py b/FIAT/arnold_winther.py index 5c12fcc76..e7af81783 100644 --- a/FIAT/arnold_winther.py +++ b/FIAT/arnold_winther.py @@ -9,157 +9,125 @@ # 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, polynomial_set +from FIAT.reference_element import TRIANGLE +from FIAT.quadrature_schemes import create_quadrature +from FIAT.functional import (ComponentPointEvaluation, + TensorBidirectionalIntegralMoment, + 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() + 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, 6), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] - dof_ids[1][entity_id] = list(range(dof_cur, dof_cur+4)) - dof_cur += 4 + # 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, 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 - 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) + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) + Q = create_quadrature(ref_el, degree) + 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))) # 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, qdegree)) + 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, formdegree, mapping=mapping) - super().__init__(Ps, Ls, degree, 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 += [IntegralLegendreNormalNormalMoment(cell, entity_id, order, 6), - IntegralLegendreNormalTangentialMoment(cell, entity_id, order, 6)] - 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))) - 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} + max_order = degree - 2 + qdegree = degree + max_order + for entity in sorted(top[1]): + cur = len(nodes) + 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} + 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, scale="L2 piola") + phis = P.tabulate(Q.get_points())[(0,)*sd] + 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} + 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) + + 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/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 1284cd5e2..6a831a87c 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): @@ -169,7 +167,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())]} super().__init__(ref_el, tuple(), pt_dict, {}, "PointEval") def __call__(self, fn): @@ -193,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): @@ -296,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): @@ -338,23 +336,41 @@ 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, 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]: + 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, {}, 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=""): - 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)} + 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, quad_deg, entity=entity) + x = cell.compute_barycentric_coordinates(Q.get_points(), entity=entity) - super().__init__(cell, shp, pt_dict, {}, nm) + 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, nm=nm) class IntegralLegendreNormalMoment(IntegralLegendreDirectionalMoment): @@ -373,38 +389,17 @@ 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): """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") @@ -412,12 +407,20 @@ 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") +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_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 against some tabulated function f.""" @@ -462,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.""" @@ -581,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): @@ -658,33 +620,25 @@ 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") -class TensorBidirectionalMomentInnerProductEvaluation(Functional): +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 \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, 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") + def __init__(self, ref_el, v, w, Q, 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") class IntegralMomentOfNormalEvaluation(Functional): @@ -730,27 +684,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") 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/hellan_herrmann_johnson.py b/FIAT/hellan_herrmann_johnson.py index 631272c64..d90a9c19f 100644 --- a/FIAT/hellan_herrmann_johnson.py +++ b/FIAT/hellan_herrmann_johnson.py @@ -3,79 +3,100 @@ # Copyright (C) 2016-2018 Lizao Li # +# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2024 +# # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later - from FIAT import dual_set, finite_element, polynomial_set -from FIAT.functional import ComponentPointEvaluation, PointwiseInnerProductEvaluation as InnerProduct +from FIAT.check_format_variant import check_format_variant +from FIAT.functional import (PointwiseInnerProductEvaluation, + ComponentPointEvaluation, + TensorBidirectionalIntegralMoment as BidirectionalMoment) +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature class HellanHerrmannJohnsonDual(dual_set.DualSet): - """Degrees of freedom for Hellan-Herrmann-Johnson elements.""" - def __init__(self, ref_el, degree): + def __init__(self, ref_el, degree, variant, qdegree): sd = ref_el.get_spatial_dimension() - if sd != 2: - raise ValueError("Hellan-Herrmann-Johnson elements are only" - "defined in dimension 2.") - top = ref_el.get_topology() + n = list(map(ref_el.compute_scaled_normal, sorted(top[sd-1]))) entity_ids = {dim: {i: [] for i in sorted(top[dim])} for dim in sorted(top)} nodes = [] - # no vertex dofs - edge_dofs, entity_ids[1] = self._generate_edge_dofs(ref_el, degree, 0) - nodes.extend(edge_dofs) - cell_nodes, entity_ids[sd] = self._generate_cell_dofs(ref_el, degree, len(nodes)) - nodes.extend(cell_nodes) + # Face dofs + if variant == "point": + # n^T u n evaluated on a Pk lattice + for f in sorted(top[sd-1]): + cur = len(nodes) + pts = ref_el.make_points(sd-1, f, degree + sd) + nodes.extend(PointwiseInnerProductEvaluation(ref_el, n[f], n[f], pt) + for pt in pts) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) - super().__init__(nodes, ref_el, entity_ids) + elif variant == "integral": + # n^T u n integrated against a basis for Pk + facet = ref_el.construct_subelement(sd-1) + Q = create_quadrature(facet, qdegree + degree) + P = polynomial_set.ONPolynomialSet(facet, degree) + phis = P.tabulate(Q.get_points())[(0,)*(sd-1)] + for f in sorted(top[sd-1]): + cur = len(nodes) + Q_mapped = FacetQuadratureRule(ref_el, sd-1, f, Q) + detJ = Q_mapped.jacobian_determinant() + nodes.extend(BidirectionalMoment(ref_el, n[f], n[f]/detJ, Q_mapped, phi) for phi in phis) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) - @staticmethod - def _generate_edge_dofs(ref_el, degree, offset): - """Generate dofs on edges. - On each edge, let n be its normal. For degree=r, the scalar function - n^T u n - is evaluated at enough points to control P(r). - """ - dim = 1 - top = ref_el.get_topology() - nodes = [] - entity_ids = {} - for entity_id in sorted(top[dim]): - pts = ref_el.make_points(dim, entity_id, degree + 2) - normal = ref_el.compute_scaled_normal(entity_id) - nodes.extend(InnerProduct(ref_el, normal, normal, pt) for pt in pts) - num_new_nodes = len(pts) - entity_ids[entity_id] = list(range(offset, offset + num_new_nodes)) - offset += num_new_nodes - return nodes, entity_ids + # Interior dofs + cur = len(nodes) + if sd == 2 and variant == "point": + # FIXME Keeping Cartesian dofs in 2D just to make regression test pass + pts = ref_el.make_points(sd, 0, degree + sd) + nodes.extend(ComponentPointEvaluation(ref_el, (i, j), (sd, sd), pt) + for i in range(sd) for j in range(i, sd) for pt in pts) + elif variant == "point": + # n[f]^T u n[f] evaluated on a P_{k-1} lattice + pts = ref_el.make_points(sd, 0, degree + sd) + nodes.extend(PointwiseInnerProductEvaluation(ref_el, n[f], n[f], pt) + for pt in pts for f in sorted(top[sd-1])) - @staticmethod - def _generate_cell_dofs(ref_el, degree, offset): - """Generate dofs on the cell interior. - On each triangle, for degree=r, the three components - u11, u12, u22 - are evaluated at enough points to control P(r-1). - """ - sd = ref_el.get_spatial_dimension() - shp = (sd, sd) - basis = [(i, j) for i in range(sd) for j in range(i, sd)] - pts = ref_el.make_points(sd, 0, degree + 2) - nodes = [ComponentPointEvaluation(ref_el, comp, shp, pt) - for comp in basis for pt in pts] - entity_ids = {0: list(range(offset, offset + len(nodes)))} - return nodes, entity_ids + # n[i+1]^T u n[i+2] evaluated on a Pk lattice + pts = ref_el.make_points(sd, 0, degree + sd + 1) + nodes.extend(PointwiseInnerProductEvaluation(ref_el, n[i+1], n[i+2], pt) + for pt in pts for i in range((sd-1)*(sd-2))) + else: + Q = create_quadrature(ref_el, qdegree + degree) + P = polynomial_set.ONPolynomialSet(ref_el, degree) + phis = P.tabulate(Q.get_points())[(0,)*sd] + phis /= ref_el.volume() + dimPkm1 = P.expansion_set.get_num_members(degree-1) + # n[f]^T u n[f] integrated against a basis for P_{k-1} + nodes.extend(BidirectionalMoment(ref_el, n[f], n[f], Q, phi) + for phi in phis[:dimPkm1] for f in sorted(top[sd-1])) + + # n[i+1]^T u n[i+2] integrated against a basis for Pk + nodes.extend(BidirectionalMoment(ref_el, n[i+1], n[i+2], Q, phi) + for phi in phis for i in range((sd-1)*(sd-2))) + + entity_ids[sd][0].extend(range(cur, len(nodes))) + + super().__init__(nodes, ref_el, entity_ids) class HellanHerrmannJohnson(finite_element.CiarletElement): - """The definition of Hellan-Herrmann-Johnson element. It is defined only in - dimension 2. It consists of piecewise polynomial symmetric-matrix-valued - functions of degree r or less with normal-normal continuity. + """The definition of Hellan-Herrmann-Johnson element. + HHJ(k) is the space of symmetric-matrix-valued polynomials of degree k + or less with normal-normal continuity. """ - def __init__(self, ref_el, degree): - assert degree >= 0, "Hellan-Herrmann-Johnson starts at degree 0!" + def __init__(self, ref_el, degree=0, variant=None): + if degree < 0: + raise ValueError(f"{type(self).__name__} only defined for degree >= 0") + + variant, qdegree = check_format_variant(variant, degree) poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) - dual = HellanHerrmannJohnsonDual(ref_el, degree) + dual = HellanHerrmannJohnsonDual(ref_el, degree, variant, qdegree) + sd = ref_el.get_spatial_dimension() + formdegree = (sd-1, sd-1) mapping = "double contravariant piola" - super().__init__(poly_set, dual, degree, mapping=mapping) + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/FIAT/hu_zhang.py b/FIAT/hu_zhang.py new file mode 100644 index 000000000..5068a6b39 --- /dev/null +++ b/FIAT/hu_zhang.py @@ -0,0 +1,89 @@ +# -*- 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 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_schemes import create_quadrature +from FIAT.functional import (ComponentPointEvaluation, + PointwiseInnerProductEvaluation, + TensorBidirectionalIntegralMoment, + IntegralLegendreNormalNormalMoment, + IntegralLegendreNormalTangentialMoment) + + +class HuZhangDual(dual_set.DualSet): + def __init__(self, ref_el, degree, variant, qdegree): + 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 + 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 + for entity in sorted(top[1]): + cur = len(nodes) + 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 + cur = len(nodes) + if variant == "point": + # 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": + # Moments of unique components against a basis for P_{k-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(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))) + 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=None): + if degree < 3: + raise ValueError(f"{type(self).__name__} only defined for degree >= 3") + if ref_el.shape != TRIANGLE: + 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) + formdegree = ref_el.get_spatial_dimension() - 1 + 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..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,33 +18,30 @@ 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) - 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 + 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, scale="L2 piola") + phis = P.tabulate(Q.get_points())[(0,) * sd] + 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))) @@ -58,5 +55,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) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 0dbc76262..2c623c64f 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -8,154 +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 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) +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) - pts = Q.get_points() - onp = onp.tabulate(pts, 0)[0, 0] + 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] + for phi in phis: + yield IntegralMomentOfDivergence(ref_el, Q, phi) - 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, :])) +class MardalTaiWintherDual(dual_set.DualSet): + """Degrees of freedom for Mardal-Tai-Winther elements.""" + def __init__(self, ref_el, degree): + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() - return ells + 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.") -class MardalTaiWintherDual(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] = dof_ids[1][entity_id] + _edge_dof_ids[entity_id] - - dof_ids[2][0] = dof_ids[2][0] + _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} + 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 = make_quadrature(facet, 6) - Pq = ONPolynomialSet(facet, 1) - Pq_at_qpts = Pq.tabulate(Q.get_points())[tuple([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)) - - 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 entity_id in range(3): - dofs += [IntegralLegendreNormalMoment(cell, entity_id, 2, 6), - IntegralLegendreNormalMoment(cell, entity_id, 3, 6)] - - edge_dof_ids[entity_id] = [offset, offset+1] - offset += 2 - - cell_dofs = DivergenceDubinerMoments(cell, 1, 2, 6) - dofs.extend(cell_dofs) - cell_dof_ids = list(range(offset, offset+len(cell_dofs))) - - return (dofs, edge_dof_ids, cell_dof_ids) - - -class MardalTaiWinther(CiarletElement): + Q = create_quadrature(facet, degree+1) + Pq = polynomial_set.ONPolynomialSet(facet, 1) + phis = Pq.tabulate(Q.get_points())[(0,)*(sd - 1)] + 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 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(nodes) + nodes.extend(DivergenceDubinerMoments(ref_el, start_order-1, stop_order-1, degree)) + entity_ids[sd][0].extend(range(cur, len(nodes))) + + 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 = 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) diff --git a/FIAT/polynomial_set.py b/FIAT/polynomial_set.py index 3a3e9ce3a..08e62b8d8 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) @@ -206,19 +203,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/FIAT/quadrature_schemes.py b/FIAT/quadrature_schemes.py index 4d48a09fe..d0c03ade6 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 @@ -27,32 +32,41 @@ # First added: 2011-04-19 # Last changed: 2011-04-19 -# NumPy import numpy -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"): +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 a + 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" -> 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 + integration domain. If not provided, the domain is the entire cell. """ + if entity is not None: + dimension, entity_id = entity + sub_el = ref_el.construct_subelement(dimension) + Q_ref = create_quadrature(sub_el, degree, scheme=scheme) + return FacetQuadratureRule(ref_el, dimension, entity_id, Q_ref) + if ref_el.is_macrocell(): dimension = ref_el.get_dimension() sub_el = ref_el.construct_subelement(dimension) diff --git a/FIAT/reference_element.py b/FIAT/reference_element.py index 822f48de1..7d9ff2968 100644 --- a/FIAT/reference_element.py +++ b/FIAT/reference_element.py @@ -626,7 +626,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/FIAT/regge.py b/FIAT/regge.py index 6e08d0399..31dd6bb44 100644 --- a/FIAT/regge.py +++ b/FIAT/regge.py @@ -3,51 +3,71 @@ # Copyright (C) 2015-2018 Lizao Li # +# Modified by Pablo D. Brubeck (brubeck@protonmail.com), 2024 +# # This file is part of FIAT (https://www.fenicsproject.org) # # SPDX-License-Identifier: LGPL-3.0-or-later from FIAT import dual_set, finite_element, polynomial_set -from FIAT.functional import PointwiseInnerProductEvaluation as InnerProduct +from FIAT.check_format_variant import check_format_variant +from FIAT.functional import (PointwiseInnerProductEvaluation, + TensorBidirectionalIntegralMoment as BidirectionalMoment) +from FIAT.quadrature import FacetQuadratureRule +from FIAT.quadrature_schemes import create_quadrature class ReggeDual(dual_set.DualSet): - """Degrees of freedom for generalized Regge finite elements. - - On a k-face for degree r, the dofs are given by the value of - t^T u t - evaluated at enough points to control P(r-k+1) for all the edge - tangents of the face. - `ref_el.make_points(dim, entity, degree + 2)` happens to - generate exactly those points needed. - """ - def __init__(self, ref_el, degree): + def __init__(self, ref_el, degree, variant, qdegree): top = ref_el.get_topology() entity_ids = {dim: {i: [] for i in sorted(top[dim])} for dim in sorted(top)} nodes = [] - for dim in sorted(top): - if dim == 0: - # no vertex dofs - continue - for entity in sorted(top[dim]): - cur = len(nodes) - tangents = ref_el.compute_face_edge_tangents(dim, entity) - pts = ref_el.make_points(dim, entity, degree + 2) - nodes.extend(InnerProduct(ref_el, t, t, pt) - for pt in pts - for t in tangents) - entity_ids[dim][entity].extend(range(cur, len(nodes))) + if variant == "point": + # On a dim-facet, for all the edge tangents of the facet, + # t^T u t is evaluated on a Pk lattice, where k = degree - dim + 1. + for dim in sorted(top): + for entity in sorted(top[dim]): + cur = len(nodes) + tangents = ref_el.compute_face_edge_tangents(dim, entity) + pts = ref_el.make_points(dim, entity, degree + 2) + nodes.extend(PointwiseInnerProductEvaluation(ref_el, t, t, pt) + for pt in pts for t in tangents) + entity_ids[dim][entity].extend(range(cur, len(nodes))) + + elif variant == "integral": + # On a dim-facet, for all the edge tangents of the facet, + # t^T u t is integrated against a basis for Pk, where k = degree - dim + 1. + for dim in sorted(top): + k = degree - dim + 1 + if dim == 0 or k < 0: + continue + facet = ref_el.construct_subelement(dim) + Q = create_quadrature(facet, qdegree + k) + P = polynomial_set.ONPolynomialSet(facet, k) + phis = P.tabulate(Q.get_points())[(0,)*dim] + for entity in sorted(top[dim]): + cur = len(nodes) + tangents = ref_el.compute_face_edge_tangents(dim, entity) + Q_mapped = FacetQuadratureRule(ref_el, dim, entity, Q) + detJ = Q_mapped.jacobian_determinant() + nodes.extend(BidirectionalMoment(ref_el, t, t/detJ, Q_mapped, phi) + for phi in phis for t in tangents) + entity_ids[dim][entity].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) class Regge(finite_element.CiarletElement): """The generalized Regge elements for symmetric-matrix-valued functions. - REG(r) in dimension n is the space of polynomial symmetric-matrix-valued - functions of degree r or less with tangential-tangential continuity. + REG(k) is the space of symmetric-matrix-valued polynomials of degree k + or less with tangential-tangential continuity. """ - def __init__(self, ref_el, degree): - assert degree >= 0, "Regge start at degree 0!" + def __init__(self, ref_el, degree=0, variant=None): + if degree < 0: + raise ValueError(f"{type(self).__name__} only defined for degree >= 0") + + variant, qdegree = check_format_variant(variant, degree) poly_set = polynomial_set.ONSymTensorPolynomialSet(ref_el, degree) - dual = ReggeDual(ref_el, degree) + dual = ReggeDual(ref_el, degree, variant, qdegree) + formdegree = (1, 1) mapping = "double covariant piola" - super().__init__(poly_set, dual, degree, mapping=mapping) + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/test/regression/test_regression.py b/test/regression/test_regression.py index af91ebded..a4101d0cb 100644 --- a/test/regression/test_regression.py +++ b/test/regression/test_regression.py @@ -280,11 +280,14 @@ def test_quadrature(quadrature_reference_data, quadrature_test_case): def create_data(family, dim, degree): '''Create the reference data. ''' + kwargs = {} + if family in {"Regge", "Hellan-Herrmann-Johnson"}: + kwargs["variant"] = "point" # Get domain and element class domain = ufc_simplex(dim) ElementClass = supported_elements[family] # Create element - element = ElementClass(domain, degree) + element = ElementClass(domain, degree, **kwargs) # Create quadrature points quad_rule = make_quadrature(domain, num_points) points = quad_rule.get_points() diff --git a/test/unit/test_awc.py b/test/unit/test_awc.py index f40690c48..b96bc86be 100644 --- a/test/unit/test_awc.py +++ b/test/unit/test_awc.py @@ -1,11 +1,11 @@ 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(): 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 @@ -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): @@ -90,11 +96,11 @@ 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) - 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])) diff --git a/test/unit/test_fiat.py b/test/unit/test_fiat.py index fb9122cb8..34d5e333d 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 @@ -42,6 +44,9 @@ 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 from FIAT.argyris import Argyris # noqa: F401 from FIAT.hermite import CubicHermite # noqa: F401 @@ -261,9 +266,28 @@ def __init__(self, a, b): "Regge(S, 0)", "Regge(S, 1)", "Regge(S, 2)", + "Regge(T, 1, variant='point')", + "Regge(S, 1, variant='point')", "HellanHerrmannJohnson(T, 0)", "HellanHerrmannJohnson(T, 1)", "HellanHerrmannJohnson(T, 2)", + "HellanHerrmannJohnson(S, 0)", + "HellanHerrmannJohnson(S, 1)", + "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)", @@ -309,6 +333,13 @@ def __init__(self, a, b): "Morley(T)", "BernardiRaugel(T)", "BernardiRaugel(S)", + "MardalTaiWinther(T, 3)", + "ArnoldWintherNC(T, 2)", + "ArnoldWinther(T, 3)", + "HuZhang(T, 3)", + "HuZhang(T, 4)", + "HuZhang(T, 3, 'point')", + "HuZhang(T, 4, 'point')", "KongMulderVeldhuizen(T,1)", "KongMulderVeldhuizen(T,2)", "KongMulderVeldhuizen(T,3)", @@ -582,98 +613,6 @@ def test_error_point_high_order(element): eval(element) -@pytest.mark.parametrize('cell', [I, T, S]) -def test_expansion_orthonormality(cell): - from FIAT import expansions - from FIAT.quadrature_schemes import create_quadrature - U = expansions.ExpansionSet(cell) - degree = 10 - rule = create_quadrature(cell, 2*degree) - phi = U.tabulate(degree, rule.pts) - qwts = rule.get_weights() - scale = 2 ** cell.get_spatial_dimension() - results = scale * np.dot(np.multiply(phi, qwts), phi.T) - - assert np.allclose(results, np.diag(np.diag(results))) - assert np.allclose(np.diag(results), 1.0) - - -@pytest.mark.parametrize('dim', range(1, 4)) -def test_expansion_values(dim): - import sympy - from FIAT import expansions, polynomial_set, reference_element - cell = reference_element.default_simplex(dim) - U = expansions.ExpansionSet(cell) - dpoints = [] - rpoints = [] - - npoints = 4 - interior = 1 - for alpha in reference_element.lattice_iter(interior, npoints+1-interior, dim): - dpoints.append(tuple(2*np.array(alpha, dtype="d")/npoints-1)) - rpoints.append(tuple(2*sympy.Rational(a, npoints)-1 for a in alpha)) - - n = 16 - Uvals = U.tabulate(n, dpoints) - idx = (lambda p: p, expansions.morton_index2, expansions.morton_index3)[dim-1] - eta = sympy.DeferredVector("eta") - half = sympy.Rational(1, 2) - - def duffy_coords(pt): - if len(pt) == 1: - return pt - elif len(pt) == 2: - eta0 = 2 * (1 + pt[0]) / (1 - pt[1]) - 1 - eta1 = pt[1] - return eta0, eta1 - else: - eta0 = 2 * (1 + pt[0]) / (-pt[1] - pt[2]) - 1 - eta1 = 2 * (1 + pt[1]) / (1 - pt[2]) - 1 - eta2 = pt[2] - return eta0, eta1, eta2 - - def basis(dim, p, q=0, r=0): - if dim >= 1: - f = sympy.jacobi(p, 0, 0, eta[0]) - f *= sympy.sqrt(half + p) - if dim >= 2: - f *= sympy.jacobi(q, 2*p+1, 0, eta[1]) * ((1 - eta[1])/2) ** p - f *= sympy.sqrt(1 + p + q) - if dim >= 3: - f *= sympy.jacobi(r, 2*p+2*q+2, 0, eta[2]) * ((1 - eta[2])/2) ** (p+q) - f *= sympy.sqrt(1 + half + p + q + r) - return f - - def eval_basis(f, pt): - return float(f.subs(dict(zip(eta, duffy_coords(pt))))) - - for i in range(n + 1): - for indices in polynomial_set.mis(dim, i): - phi = basis(dim, *indices) - exact = np.array([eval_basis(phi, r) for r in rpoints]) - uh = Uvals[idx(*indices)] - assert np.allclose(uh, exact, atol=1E-14) - - -@pytest.mark.parametrize('cell', [I, T, S]) -def test_bubble_duality(cell): - from FIAT.polynomial_set import make_bubbles - from FIAT.quadrature_schemes import create_quadrature - degree = 10 - sd = cell.get_spatial_dimension() - B = make_bubbles(cell, degree) - - Q = create_quadrature(cell, 2*B.degree - sd - 1) - qpts, qwts = Q.get_points(), Q.get_weights() - phi = B.tabulate(qpts)[(0,) * sd] - phi_dual = phi / abs(phi[0]) - scale = 2 ** sd - results = scale * np.dot(np.multiply(phi_dual, qwts), phi.T) - - assert np.allclose(results, np.diag(np.diag(results))) - assert np.allclose(np.diag(results), 1.0) - - if __name__ == '__main__': import os pytest.main(os.path.abspath(__file__)) 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) diff --git a/test/unit/test_polynomial.py b/test/unit/test_polynomial.py new file mode 100644 index 000000000..2982cf11e --- /dev/null +++ b/test/unit/test_polynomial.py @@ -0,0 +1,109 @@ +# Copyright (C) 2024 Pablo Brubeck +# +# This file is part of FIAT. +# +# FIAT is free software: you can redistribute it and/or modify +# it under the terms of the GNU Lesser General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# FIAT is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU Lesser General Public License for more details. +# +# You should have received a copy of the GNU Lesser General Public License +# along with FIAT. If not, see . + +import pytest +import numpy +import sympy + +from FIAT import expansions, polynomial_set, reference_element +from FIAT.quadrature_schemes import create_quadrature + + +@pytest.fixture(params=(1, 2, 3)) +def cell(request): + dim = request.param + return reference_element.default_simplex(dim) + + +@pytest.mark.parametrize("degree", [10]) +def test_expansion_values(cell, degree): + dim = cell.get_spatial_dimension() + U = expansions.ExpansionSet(cell) + dpoints = [] + rpoints = [] + + numpyoints = 4 + interior = 1 + for alpha in reference_element.lattice_iter(interior, numpyoints+1-interior, dim): + dpoints.append(tuple(2*numpy.array(alpha, dtype="d")/numpyoints-1)) + rpoints.append(tuple(2*sympy.Rational(a, numpyoints)-1 for a in alpha)) + + Uvals = U.tabulate(degree, dpoints) + idx = (lambda p: p, expansions.morton_index2, expansions.morton_index3)[dim-1] + eta = sympy.DeferredVector("eta") + half = sympy.Rational(1, 2) + + def duffy_coords(pt): + if len(pt) == 1: + return pt + elif len(pt) == 2: + eta0 = 2 * (1 + pt[0]) / (1 - pt[1]) - 1 + eta1 = pt[1] + return eta0, eta1 + else: + eta0 = 2 * (1 + pt[0]) / (-pt[1] - pt[2]) - 1 + eta1 = 2 * (1 + pt[1]) / (1 - pt[2]) - 1 + eta2 = pt[2] + return eta0, eta1, eta2 + + def basis(dim, p, q=0, r=0): + if dim >= 1: + f = sympy.jacobi(p, 0, 0, eta[0]) + f *= sympy.sqrt(half + p) + if dim >= 2: + f *= sympy.jacobi(q, 2*p+1, 0, eta[1]) * ((1 - eta[1])/2) ** p + f *= sympy.sqrt(1 + p + q) + if dim >= 3: + f *= sympy.jacobi(r, 2*p+2*q+2, 0, eta[2]) * ((1 - eta[2])/2) ** (p+q) + f *= sympy.sqrt(1 + half + p + q + r) + return f + + def eval_basis(f, pt): + return float(f.subs(dict(zip(eta, duffy_coords(pt))))) + + for i in range(degree + 1): + for indices in polynomial_set.mis(dim, i): + phi = basis(dim, *indices) + exact = numpy.array([eval_basis(phi, r) for r in rpoints]) + uh = Uvals[idx(*indices)] + assert numpy.allclose(uh, exact, atol=1E-14) + + +@pytest.mark.parametrize("degree", [10]) +def test_expansion_orthonormality(cell, degree): + U = expansions.ExpansionSet(cell) + rule = create_quadrature(cell, 2*degree) + phi = U.tabulate(degree, rule.pts) + qwts = rule.get_weights() + results = numpy.dot(numpy.multiply(phi, qwts), phi.T) + assert numpy.allclose(results, numpy.diag(numpy.diag(results))) + assert numpy.allclose(numpy.diag(results), 1.0) + + +@pytest.mark.parametrize("degree", [10]) +def test_bubble_duality(cell, degree): + sd = cell.get_spatial_dimension() + B = polynomial_set.make_bubbles(cell, degree) + + Q = create_quadrature(cell, 2*B.degree - sd - 1) + qpts, qwts = Q.get_points(), Q.get_weights() + phi = B.tabulate(qpts)[(0,) * sd] + phi_dual = phi / abs(phi[0]) + scale = 2 ** sd + results = scale * numpy.dot(numpy.multiply(phi_dual, qwts), phi.T) + assert numpy.allclose(results, numpy.diag(numpy.diag(results))) + assert numpy.allclose(numpy.diag(results), 1.0)