diff --git a/FIAT/__init__.py b/FIAT/__init__.py index 4b90bcb2..5ecfba55 100644 --- a/FIAT/__init__.py +++ b/FIAT/__init__.py @@ -27,7 +27,7 @@ from FIAT.arnold_qin import ArnoldQin from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanH1div from FIAT.christiansen_hu import ChristiansenHu -from FIAT.johnson_mercier import JohnsonMercier +from FIAT.johnson_mercier import JohnsonMercier, ReducedJohnsonMercier from FIAT.brezzi_douglas_marini import BrezziDouglasMarini from FIAT.Sminus import TrimmedSerendipityEdge, TrimmedSerendipityFace from FIAT.SminusDiv import TrimmedSerendipityDiv @@ -104,6 +104,7 @@ "Guzman-Neilan 2nd kind H1": GuzmanNeilanSecondKindH1, "Guzman-Neilan H1(div)": GuzmanNeilanH1div, "Johnson-Mercier": JohnsonMercier, + "Reduced Johnson-Mercier": ReducedJohnsonMercier, "Lagrange": Lagrange, "Kong-Mulder-Veldhuizen": KongMulderVeldhuizen, "Gauss-Lobatto-Legendre": GaussLobattoLegendre, diff --git a/FIAT/bernardi_raugel.py b/FIAT/bernardi_raugel.py index 852f2b0d..ea2fac3a 100644 --- a/FIAT/bernardi_raugel.py +++ b/FIAT/bernardi_raugel.py @@ -116,6 +116,9 @@ def __init__(self, ref_el, order=1, degree=None, reduced=False, ref_complex=None Q = Qt[f] phi = ft_at_qpts udir = thats[f][i-1] + if sd == 3: + udir = numpy.cross(perp(*thats[f]), udir) + nodes.append(FrobeniusIntegralMoment(ref_el, Q, numpy.outer(udir, phi))) entity_ids[sd-1][f].extend(range(cur, len(nodes))) super().__init__(nodes, ref_el, entity_ids) diff --git a/FIAT/functional.py b/FIAT/functional.py index 82509058..4041b73e 100644 --- a/FIAT/functional.py +++ b/FIAT/functional.py @@ -459,6 +459,7 @@ def __init__(self, ref_el, Q, f_at_qpts): self.Q = Q sd = ref_el.get_spatial_dimension() + shp = f_at_qpts.shape[1:-1] + (sd,) points = Q.get_points() self.dpts = points @@ -468,7 +469,7 @@ def __init__(self, ref_el, Q, f_at_qpts): dpt_dict = {tuple(pt): [(wt, alphas[i], (i,)) for i in range(sd)] for pt, wt in zip(points, weights)} - super().__init__(ref_el, tuple(), {}, dpt_dict, + super().__init__(ref_el, shp, {}, dpt_dict, "IntegralMomentOfDivergence") @@ -492,7 +493,7 @@ def __init__(self, ref_el, Q, f_at_qpts): dpt_dict = {tuple(pt): [(wt[i], alphas[j], (i, j)) for i, j in numpy.ndindex(shp)] for pt, wt in zip(points, weights)} - super().__init__(ref_el, tuple(), {}, dpt_dict, "IntegralMomentOfDivergence") + super().__init__(ref_el, shp, {}, dpt_dict, "IntegralMomentOfDivergence") class PointNormalEvaluation(Functional): diff --git a/FIAT/johnson_mercier.py b/FIAT/johnson_mercier.py index 4521c016..2b20fcf4 100644 --- a/FIAT/johnson_mercier.py +++ b/FIAT/johnson_mercier.py @@ -1,7 +1,9 @@ from FIAT import finite_element, dual_set, macro, polynomial_set from FIAT.check_format_variant import parse_quadrature_scheme -from FIAT.functional import TensorBidirectionalIntegralMoment +from FIAT.functional import (FrobeniusIntegralMoment, IntegralMomentOfTensorDivergence, + TensorBidirectionalIntegralMoment) from FIAT.quadrature import FacetQuadratureRule +from FIAT.nedelec_second_kind import NedelecSecondKind import numpy @@ -28,7 +30,12 @@ def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): cur = len(nodes) Q = FacetQuadratureRule(ref_el, dim, f, Qref, avg=True) thats = ref_el.compute_tangents(dim, f) - nhat = numpy.dot(R, *thats) if sd == 2 else numpy.cross(*thats) + if sd == 2: + nhat = numpy.dot(R, *thats) + else: + nhat = numpy.cross(*thats) + thats = numpy.cross(nhat[None, :], thats, axis=1) + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, nhat, comp, Q, phi) for phi in phis for comp in (nhat, *thats)) entity_ids[dim][f].extend(range(cur, len(nodes))) @@ -57,3 +64,124 @@ def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): formdegree = ref_el.get_spatial_dimension() - 1 mapping = "double contravariant piola" super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) + + +class DG0Alfeld(finite_element.CiarletElement): + + def __init__(self, ref_complex, P1=None, quad_scheme=None): + ref_el = ref_complex.get_parent() + sd = ref_el.get_spatial_dimension() + if P1 is None: + if sd == 1: + P1 = polynomial_set.ONPolynomialSet(ref_el, 1, shape=(sd,)) + else: + P1 = NedelecSecondKind(ref_el, 1) + + Q = parse_quadrature_scheme(ref_complex, P1.degree(), quad_scheme) + phis = P1.tabulate(0, Q.get_points())[(0,)*sd] + nodes = [FrobeniusIntegralMoment(ref_el, Q, phi) for phi in phis] + dual = dual_set.DualSet(nodes, ref_el, P1.entity_dofs()) + + poly_set = polynomial_set.ONPolynomialSet(ref_complex, 0, shape=(sd,)) + formdegree = sd + super().__init__(poly_set, dual, poly_set.degree, formdegree) + + +def rbm_complement(ref_complex): + """Constructs a basis for the complement of the rigid body motions over P0(Alfeld).""" + P0 = DG0Alfeld(ref_complex) + P0 = DG0Alfeld(ref_complex, P1=P0) + + entity_ids = P0.entity_dofs() + ids = [] + for entity in entity_ids[1]: + ids.extend(entity_ids[1][entity][1:]) + return P0.get_nodal_basis().take(ids) + + +class ReducedJohnsonMercierDualSet(dual_set.DualSet): + def __init__(self, ref_complex, degree, variant=None, quad_scheme=None): + if degree != 1: + raise ValueError("Reduced Johnson-Mercier only defined for degree=1") + if variant is not None: + raise ValueError(f"Reduced Johnson-Mercier does not have the {variant} variant") + ref_el = ref_complex.get_parent() + sd = ref_el.get_spatial_dimension() + top = ref_el.get_topology() + entity_ids = {dim: {entity: [] for entity in sorted(top[dim])} for dim in sorted(top)} + nodes = [] + + ref_facet = ref_el.get_facet_element() + Q = parse_quadrature_scheme(ref_facet, degree+1, quad_scheme) + P1 = polynomial_set.ONPolynomialSet(ref_facet, degree, scale="orthonormal") + P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] + dimP1 = len(P1)*(sd-1) + dimNed1 = dimP1 // 2 + if sd == 3: + # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] + RT_at_qpts = numpy.zeros((dimP1, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0] = P1_at_qpts[0] + RT_at_qpts[1, 1] = P1_at_qpts[0] + RT_at_qpts[2, 0] = P1_at_qpts[1] + RT_at_qpts[2, 1] = P1_at_qpts[2] + # Basis for the complement of RT [(y, x), (x, -y), (y, -x)] + RT_at_qpts[3, 0] = P1_at_qpts[2] + RT_at_qpts[3, 1] = P1_at_qpts[1] + RT_at_qpts[4, 0] = P1_at_qpts[1] + RT_at_qpts[4, 1] = -P1_at_qpts[2] + RT_at_qpts[5, 0] = P1_at_qpts[2] + RT_at_qpts[5, 1] = -P1_at_qpts[1] + else: + RT_at_qpts = P1_at_qpts[:, None, :] + + # Facet degrees of freedom + for f in sorted(top[sd-1]): + cur = len(nodes) + n = ref_el.compute_scaled_normal(f) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = Qf.jacobian() + # Normal-normal moments against P1 + nodes.extend(TensorBidirectionalIntegralMoment(ref_el, n, n, Qf, phi) for phi in P1_at_qpts) + # Normal-tangential moments against n x RT + phis = numpy.tensordot(Jf, RT_at_qpts[:dimNed1].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + phis = numpy.cross(n[None, :, None], phis, axis=1) + phis = phis[:, :, None, :] * n[None, None, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + # Facet constraints + for f in sorted(top[sd-1]): + cur = len(nodes) + n = ref_el.compute_scaled_normal(f) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + Jf = Qf.jacobian() + # Normal-tangential moments against n x (P1 \ RT) + phis = numpy.tensordot(Jf, RT_at_qpts[dimNed1:].transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + phis = numpy.cross(n[None, :, None], phis, axis=1) + phis = phis[:, :, None, :] * n[None, None, :, None] + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) + entity_ids[sd-1][f].extend(range(cur, len(nodes))) + + # Interior constraints: moments of divergence against (P1 \ Nedelec) + ref_complex = macro.AlfeldSplit(ref_el) + Q = parse_quadrature_scheme(ref_complex, 2*(degree-1)) + comp_space = rbm_complement(ref_complex) + phis = comp_space.tabulate(Q.get_points())[(0,)*sd] + cur = len(nodes) + nodes.extend(IntegralMomentOfTensorDivergence(ref_el, Q, phi) for phi in phis) + entity_ids[sd][0].extend(range(cur, len(nodes))) + super().__init__(nodes, ref_el, entity_ids) + + +class ReducedJohnsonMercier(finite_element.CiarletElement): + """The Reduced Johnson-Mercier finite element.""" + + def __init__(self, ref_el, degree=1, variant=None, quad_scheme=None): + ref_complex = macro.AlfeldSplit(ref_el) + poly_set = macro.HDivSymPolynomialSet(ref_complex, degree) + dual = ReducedJohnsonMercierDualSet(ref_complex, degree, variant=variant, quad_scheme=quad_scheme) + formdegree = ref_el.get_spatial_dimension() - 1 + mapping = "double contravariant piola" + super().__init__(poly_set, dual, degree, formdegree, mapping=mapping) diff --git a/FIAT/mardal_tai_winther.py b/FIAT/mardal_tai_winther.py index 2c623c64..7c2cf58c 100644 --- a/FIAT/mardal_tai_winther.py +++ b/FIAT/mardal_tai_winther.py @@ -8,26 +8,52 @@ # SPDX-License-Identifier: LGPL-3.0-or-later +import numpy from FIAT import dual_set, expansions, finite_element, polynomial_set -from FIAT.functional import (IntegralMomentOfNormalEvaluation, - IntegralMomentOfTangentialEvaluation, - IntegralLegendreNormalMoment, - IntegralMomentOfDivergence) - +from FIAT.functional import FrobeniusIntegralMoment +from FIAT.quadrature import FacetQuadratureRule from FIAT.quadrature_schemes import create_quadrature -def DivergenceDubinerMoments(ref_el, start_deg, stop_deg, comp_deg): +def curl(table_u): + grad_u = {alpha.index(1): table_u[alpha] for alpha in table_u if sum(alpha) == 1} + d = len(grad_u) + if d == 2: + curl_u = [grad_u[1], -grad_u[0]] + return numpy.concatenate(curl_u, axis=1) + else: + indices = ((i, j) for i in reversed(range(d)) for j in reversed(range(i+1, d))) + curl_u = [((-1)**k) * (grad_u[j][:, i, :] - grad_u[i][:, j, :]) for k, (i, j) in enumerate(indices)] + return numpy.transpose(curl_u, (1, 0, 2)) + + +def MardalTaiWintherSpace(ref_el): + """Construct the MTW space [P1]^d + curl(B [P1]^d).""" sd = ref_el.get_spatial_dimension() - P = polynomial_set.ONPolynomialSet(ref_el, stop_deg) - Q = create_quadrature(ref_el, comp_deg + stop_deg) + # Polynomials of degree sd+1 + degree = sd + 1 + Pk = polynomial_set.ONPolynomialSet(ref_el, degree, shape=(sd,), scale="orthonormal") + + # Grab [P1]^d from [Pk]^d + dimP1 = expansions.polynomial_dimension(ref_el, 1) + dimPk = expansions.polynomial_dimension(ref_el, degree) + ids = [i+dimPk*j for i in range(dimP1) for j in range(sd)] + P1 = Pk.take(ids) + # Project curl(B [P1]^d) into [Pk]^d + BP1 = polynomial_set.make_bubbles(ref_el, degree+1, shape=((sd*(sd-1))//2,)) - 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) + Q = create_quadrature(ref_el, 2*degree) + qpts = Q.get_points() + qwts = Q.get_weights() + Pk_at_qpts = Pk.tabulate(qpts) + BP1_at_qpts = BP1.tabulate(qpts, 1) + + inner = lambda u, v, qwts: numpy.tensordot(u, numpy.multiply(v, qwts), axes=(range(1, u.ndim),)*2) + C = inner(curl(BP1_at_qpts), Pk_at_qpts[(0,)*sd], qwts) + coeffs = numpy.tensordot(C, Pk.get_coeffs(), axes=(1, 0)) + curlBP1 = polynomial_set.PolynomialSet(ref_el, degree, degree, Pk.get_expansion_set(), coeffs) + + return polynomial_set.polynomial_set_union_normalized(P1, curlBP1) class MardalTaiWintherDual(dual_set.DualSet): @@ -36,58 +62,47 @@ 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 sd not in (2, 3): + raise ValueError("Mardal-Tai-Winther elements are only defined in dimension 2 and 3.") - if degree != 3: - raise ValueError("Mardal-Tai-Winther elements are only defined for degree 3.") + if degree != sd+1: + raise ValueError("Mardal-Tai-Winther elements are only defined for degree = dim+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 = 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))) + # u.n against a Dubiner basis for P1 + # and u x n against a basis for lowest-order RT. + ref_facet = ref_el.get_facet_element() + Q = create_quadrature(ref_facet, degree+1) + P1 = polynomial_set.ONPolynomialSet(ref_facet, 1) + P1_at_qpts = P1.tabulate(Q.get_points())[(0,)*(sd - 1)] + if sd == 2: + # For 2D just take the constant + RT_at_qpts = P1_at_qpts[:1, None, :] + else: + # Basis for lowest-order RT [(1, 0), (0, 1), (x, y)] + RT_at_qpts = numpy.zeros((3, sd-1, P1_at_qpts.shape[-1])) + RT_at_qpts[0, 0, :] = P1_at_qpts[0, None, :] + RT_at_qpts[1, 1, :] = P1_at_qpts[0, None, :] + RT_at_qpts[2, 0, :] = P1_at_qpts[1, None, :] + RT_at_qpts[2, 1, :] = P1_at_qpts[2, None, :] - # 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)) + n = ref_el.compute_scaled_normal(f) + Qf = FacetQuadratureRule(ref_el, sd-1, f, Q, avg=True) + # Normal moments against P1 + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, numpy.outer(n, phi)) for phi in P1_at_qpts) + # Map the RT basis into the facet + Jf = Qf.jacobian() + phis = numpy.tensordot(Jf, RT_at_qpts.transpose(1, 0, 2), (1, 0)).transpose(1, 0, 2) + if sd == 3: + # Moments against cross(n, RT) + phis = numpy.cross(n[None, :, None], phis, axis=1) + nodes.extend(FrobeniusIntegralMoment(ref_el, Qf, phi) for phi in phis) 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) @@ -95,11 +110,7 @@ class MardalTaiWinther(finite_element.CiarletElement): """The definition of the Mardal-Tai-Winther element. """ def __init__(self, ref_el, degree=3): - sd = ref_el.get_spatial_dimension() - assert degree == 3, "Only defined for degree 3" - 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().__init__(poly_set, dual, degree, formdegree, mapping=mapping) + poly_set = MardalTaiWintherSpace(ref_el) + formdegree = ref_el.get_spatial_dimension() - 1 + super().__init__(poly_set, dual, degree, formdegree, mapping="contravariant piola") diff --git a/finat/__init__.py b/finat/__init__.py index bde435a1..186816d8 100644 --- a/finat/__init__.py +++ b/finat/__init__.py @@ -28,7 +28,7 @@ from .guzman_neilan import GuzmanNeilanFirstKindH1, GuzmanNeilanSecondKindH1, GuzmanNeilanBubble, GuzmanNeilanH1div # noqa: F401 from .powell_sabin import QuadraticPowellSabin6, QuadraticPowellSabin12 # noqa: F401 from .hermite import Hermite # noqa: F401 -from .johnson_mercier import JohnsonMercier # noqa: F401 +from .johnson_mercier import JohnsonMercier, ReducedJohnsonMercier # noqa: F401 from .mtw import MardalTaiWinther # noqa: F401 from .morley import Morley # noqa: F401 from .walkington import Walkington # noqa: F401 diff --git a/finat/aw.py b/finat/aw.py index 97ad614b..f498a61b 100644 --- a/finat/aw.py +++ b/finat/aw.py @@ -28,10 +28,14 @@ def _facet_transform(fiat_cell, facet_moment_degree, coordinate_mapping): transform = normal_tangential_face_transform for f in range(num_facets): - rows = transform(fiat_cell, J, detJ, f) + *Bnt, Btt = transform(fiat_cell, J, detJ, f) for i in range(dimPk_facet): s = dofs_per_facet*f + i * sd - V[s+1:s+sd, s:s+sd] = rows + ndof = s + tdofs = range(s+1, s+sd) + V[tdofs, ndof] = Bnt + V[tdofs, tdofs] = Btt + return V diff --git a/finat/element_factory.py b/finat/element_factory.py index 8563aa60..7129200b 100644 --- a/finat/element_factory.py +++ b/finat/element_factory.py @@ -66,6 +66,7 @@ "Guzman-Neilan H1(div)": finat.GuzmanNeilanH1div, "Guzman-Neilan Bubble": finat.GuzmanNeilanBubble, "Johnson-Mercier": finat.JohnsonMercier, + "Reduced Johnson-Mercier": finat.ReducedJohnsonMercier, "Lagrange": finat.Lagrange, "Kong-Mulder-Veldhuizen": finat.KongMulderVeldhuizen, "Gauss-Lobatto-Legendre": finat.GaussLobattoLegendre, diff --git a/finat/johnson_mercier.py b/finat/johnson_mercier.py index 7b6d485c..de6c3a0d 100644 --- a/finat/johnson_mercier.py +++ b/finat/johnson_mercier.py @@ -5,25 +5,90 @@ from finat.citations import cite from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement +from finat.piola_mapped import normal_tangential_edge_transform, normal_tangential_face_transform +from copy import deepcopy class JohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued def __init__(self, cell, degree=1, variant=None, quad_scheme=None): cite("Gopalakrishnan2024") - self._indices = slice(None, None) super().__init__(FIAT.JohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) def basis_transformation(self, coordinate_mapping): - numbf = self._element.space_dimension() - ndof = self.space_dimension() - - V = identity(numbf, ndof) + V = identity(self.space_dimension()) Vsub = _facet_transform(self.cell, 1, coordinate_mapping) - Vsub = Vsub[:, self._indices] m, n = Vsub.shape V[:m, :n] = Vsub + return ListTensor(V.T) + + +class ReducedJohnsonMercier(PhysicallyMappedElement, FiatElement): # symmetric matrix valued + def __init__(self, cell, degree=1, variant=None, quad_scheme=None): + cite("Gopalakrishnan2024") + super().__init__(FIAT.ReducedJohnsonMercier(cell, degree, variant=variant, quad_scheme=quad_scheme)) + # On each facet we expect the normal dofs followed by the tangential ones + # The tangential and interior constraints should be numbered last. + sd = self.cell.get_spatial_dimension() + reduced_dofs = deepcopy(self._element.entity_dofs()) + reduced_dim = 0 + + dimP1 = sd + dimNed1 = (sd*(sd-1))//2 + r = dimP1 + dimNed1 + for entity in sorted(reduced_dofs[sd-1]): + reduced_dim += len(reduced_dofs[sd-1][entity][r:]) + reduced_dofs[sd-1][entity] = reduced_dofs[sd-1][entity][:r] + + for entity in sorted(reduced_dofs[sd]): + reduced_dim += len(reduced_dofs[sd][entity]) + reduced_dofs[sd][entity] = [] + + self._entity_dofs = reduced_dofs + self._space_dimension = self._element.space_dimension() - reduced_dim + + def entity_dofs(self): + return self._entity_dofs + + def space_dimension(self): + return self._space_dimension + + def basis_transformation(self, coordinate_mapping): + sd = self.cell.get_spatial_dimension() + bary, = self.cell.make_points(sd, 0, sd+1) + J = coordinate_mapping.jacobian_at(bary) + detJ = coordinate_mapping.detJ_at(bary) + + num_dofs = self.space_dimension() + num_bfs = self._element.space_dimension() + V = identity(num_bfs, num_dofs) + dimP1 = sd + dimNed1 = (sd*(sd-1))//2 + + if sd == 2: + transform = normal_tangential_edge_transform + else: + transform = normal_tangential_face_transform + + entity_dofs = self._element.entity_dofs() + for f in sorted(entity_dofs[sd-1]): + *Bnt, Btt = transform(self.cell, J, detJ, f) + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:(dimP1+dimNed1)] + cdofs = entity_dofs[sd-1][f][dimP1+dimNed1:] + + V[tdofs, tdofs] = Btt + if sd == 2: + V[tdofs, ndofs[0]] = Bnt + V[cdofs, ndofs[1]] = Bnt + else: + V[tdofs[:-1], ndofs[0]] = Bnt + # n x (t0 * x + t1 * y) + V[tdofs[-1], ndofs[1:]] = Bnt + # n x (t0 * y + t1 * x) + V[cdofs[0], ndofs[1:]] = (Bnt[1], Bnt[0]) + # n x (t0 * x - t1 * y) + V[cdofs[1], ndofs[1:]] = (Bnt[0], -1*Bnt[1]) + # n x (t0 * y - t1 * x) + V[cdofs[2], ndofs[1:]] = (-1*Bnt[1], Bnt[0]) - # Note: that the edge DOFs are scaled by edge lengths in FIAT implies - # that they are already have the necessary rescaling to improve - # conditioning. return ListTensor(V.T) diff --git a/finat/mtw.py b/finat/mtw.py index 9e4ebaa2..c4173e31 100644 --- a/finat/mtw.py +++ b/finat/mtw.py @@ -1,46 +1,43 @@ import FIAT - from gem import ListTensor from finat.citations import cite from finat.fiat_elements import FiatElement from finat.physically_mapped import identity, PhysicallyMappedElement -from finat.piola_mapped import normal_tangential_edge_transform -from copy import deepcopy +from finat.piola_mapped import normal_tangential_edge_transform, normal_tangential_face_transform class MardalTaiWinther(PhysicallyMappedElement, FiatElement): - def __init__(self, cell, degree=3): + def __init__(self, cell, degree=None): + if degree is None: + degree = cell.get_spatial_dimension()+1 cite("Mardal2002") super().__init__(FIAT.MardalTaiWinther(cell, degree)) - reduced_dofs = deepcopy(self._element.entity_dofs()) - sd = cell.get_spatial_dimension() - fdofs = sd + 1 - reduced_dofs[sd][0] = [] - for f in reduced_dofs[sd-1]: - reduced_dofs[sd-1][f] = reduced_dofs[sd-1][f][:fdofs] - self._entity_dofs = reduced_dofs - self._space_dimension = fdofs * len(reduced_dofs[sd-1]) - def basis_transformation(self, coordinate_mapping): - numbf = self._element.space_dimension() - ndof = self.space_dimension() - V = identity(numbf, ndof) - sd = self.cell.get_spatial_dimension() bary, = self.cell.make_points(sd, 0, sd+1) J = coordinate_mapping.jacobian_at(bary) detJ = coordinate_mapping.detJ_at(bary) + + V = identity(self.space_dimension()) + dimP1 = sd + + if sd == 2: + transform = normal_tangential_edge_transform + else: + transform = normal_tangential_face_transform + entity_dofs = self.entity_dofs() for f in sorted(entity_dofs[sd-1]): - cur = entity_dofs[sd-1][f][0] - V[cur+1, cur:cur+sd] = normal_tangential_edge_transform(self.cell, J, detJ, f) + *Bnt, Btt = transform(self.cell, J, detJ, f) + ndofs = entity_dofs[sd-1][f][:dimP1] + tdofs = entity_dofs[sd-1][f][dimP1:] + V[tdofs, tdofs] = Btt + if sd == 2: + V[tdofs, ndofs[0]] = Bnt + else: + V[tdofs[-1], ndofs[1:]] = Bnt + V[tdofs[:-1], ndofs[0]] = Bnt return ListTensor(V.T) - - def entity_dofs(self): - return self._entity_dofs - - def space_dimension(self): - return self._space_dimension diff --git a/finat/physically_mapped.py b/finat/physically_mapped.py index 978bbaef..fe004440 100644 --- a/finat/physically_mapped.py +++ b/finat/physically_mapped.py @@ -38,8 +38,8 @@ def matvec(self, table): for i, js in enumerate(self.csr)] val = gem.ListTensor(exprs) - # val = self.M @ table - return gem.optimise.aggressive_unroll(val) + # val = gem.optimise.aggressive_unroll(self.M @ table) + return val def __getitem__(self, alpha): try: diff --git a/finat/piola_mapped.py b/finat/piola_mapped.py index d5fbf8a2..40555cef 100644 --- a/finat/piola_mapped.py +++ b/finat/piola_mapped.py @@ -1,7 +1,7 @@ import numpy from finat.fiat_elements import FiatElement -from finat.physically_mapped import adjugate, identity, PhysicallyMappedElement +from finat.physically_mapped import adjugate, determinant, identity, PhysicallyMappedElement from gem import Literal, ListTensor, Zero from copy import deepcopy from itertools import chain @@ -38,21 +38,21 @@ def normal_tangential_face_transform(fiat_cell, J, detJ, f): thats = fiat_cell.compute_tangents(2, f) nhat = numpy.cross(*thats) nhat /= numpy.dot(nhat, nhat) - orth_vecs = numpy.array([nhat, - numpy.cross(nhat, thats[1]), - numpy.cross(thats[0], nhat)]) - # Compute A = (alpha, beta, gamma) - Jts = J @ Literal(thats.T) - Jorths = J @ Literal(orth_vecs.T) - A = Jorths.T @ Jts - # Compute the last two rows of inv([[1, 0, 0], A.T/detJ]) - det0 = A[1, 0] * A[2, 1] - A[1, 1] * A[2, 0] - det1 = A[2, 0] * A[0, 1] - A[2, 1] * A[0, 0] - det2 = A[0, 0] * A[1, 1] - A[0, 1] * A[1, 0] - scale = detJ / det0 - rows = ((-1 * det1 / det0, -1 * scale * A[2, 1], scale * A[2, 0]), - (-1 * det2 / det0, scale * A[1, 1], -1 * scale * A[1, 0])) - return rows + orths = numpy.cross(thats, nhat[None, :], axis=1) + + Jn = J @ Literal(nhat) + Jthats = J @ Literal(thats.T) + Jorths = J @ Literal(orths.T) + A = Jthats.T @ Jorths + B = Jn @ Jthats + A = numpy.array([[A[i, j] for j in range(A.shape[1])] for i in range(A.shape[0])]) + B = numpy.array([B[i] for i in range(B.shape[0])]) + + Q = numpy.dot(thats, thats.T) + beta = determinant(A) + alpha = Q @ (adjugate(A) @ B) + row = (alpha[0] / beta, alpha[1] / beta, detJ / beta) + return row class PiolaBubbleElement(PhysicallyMappedElement, FiatElement): @@ -118,10 +118,10 @@ def basis_transformation(self, coordinate_mapping): transform = normal_tangential_face_transform for f in sorted(dofs[sd-1]): - rows = numpy.asarray(transform(self.cell, J, detJ, f)) - cur_dofs = dofs[sd-1][f] + *Bnt, Btt = transform(self.cell, J, detJ, f) + cur_dof = dofs[sd-1][f][0] cur_bfs = bfs[sd-1][f][1:] - V[numpy.ix_(cur_bfs, cur_dofs)] = rows[..., :len(cur_dofs)] + V[cur_bfs, cur_dof] = Bnt # Fix discrepancy between normal and tangential moments needs_facet_vertex_coupling = len(dofs[0][0]) > 0 and numbf > ndof diff --git a/finat/ufl/elementlist.py b/finat/ufl/elementlist.py index 4d7ae92d..b4743094 100644 --- a/finat/ufl/elementlist.py +++ b/finat/ufl/elementlist.py @@ -111,7 +111,7 @@ def show_elements(): # Zany elements register_element("Bernardi-Raugel", "BR", 1, H1, "contravariant Piola", (1, None), simplices[1:]) register_element("Bernardi-Raugel Bubble", "BRB", 1, H1, "contravariant Piola", (None, None), simplices[1:]) -register_element("Mardal-Tai-Winther", "MTW", 1, H1, "contravariant Piola", (3, 3), ("triangle",)) +register_element("Mardal-Tai-Winther", "MTW", 1, H1, "contravariant Piola", (3, 4), ("triangle", "tetrahedron")) register_element("Hermite", "HER", 0, H1, "custom", (3, 3), simplices) register_element("Argyris", "ARG", 0, H2, "custom", (5, None), ("triangle",)) register_element("Bell", "BELL", 0, H2, "custom", (5, 5), ("triangle",)) @@ -127,6 +127,7 @@ def show_elements(): register_element("Hsieh-Clough-Tocher", "HCT", 0, H2, "custom", (3, None), ("triangle",)) register_element("Reduced-Hsieh-Clough-Tocher", "HCT-red", 0, H2, "custom", (3, 3), ("triangle",)) register_element("Johnson-Mercier", "JM", 2, HDiv, "double contravariant Piola", (1, 1), simplices[1:]) +register_element("Reduced Johnson-Mercier", "JM-red", 2, HDiv, "double contravariant Piola", (1, 1), simplices[1:]) register_element("Walkington", "WALK", 0, H2, "custom", (5, 5), ("tetrahedron",)) register_element("Alfeld C2", "ALF-C2", 0, H3, "custom", (5, None), ("triangle",)) diff --git a/test/FIAT/unit/test_fiat.py b/test/FIAT/unit/test_fiat.py index 361cd510..11f44313 100644 --- a/test/FIAT/unit/test_fiat.py +++ b/test/FIAT/unit/test_fiat.py @@ -59,7 +59,7 @@ from FIAT.christiansen_hu import ChristiansenHu # noqa: F401 from FIAT.guzman_neilan import GuzmanNeilanFirstKindH1 # noqa: F401 from FIAT.guzman_neilan import GuzmanNeilanSecondKindH1 # noqa: F401 -from FIAT.johnson_mercier import JohnsonMercier # noqa: F401 +from FIAT.johnson_mercier import JohnsonMercier, ReducedJohnsonMercier # noqa: F401 from FIAT.bubble import Bubble, FacetBubble # noqa: F401 from FIAT.enriched import EnrichedElement # noqa: F401 from FIAT.nodal_enriched import NodalEnrichedElement @@ -361,6 +361,7 @@ def __init__(self, a, b): "BernardiRaugel(T)", "BernardiRaugel(S)", "MardalTaiWinther(T, 3)", + "MardalTaiWinther(S, 4)", "ArnoldWintherNC(T, 2)", "ArnoldWinther(T, 3)", "HuZhang(T, 3)", @@ -386,6 +387,8 @@ def __init__(self, a, b): "AlfeldC2(T, 6)", "JohnsonMercier(T)", "JohnsonMercier(S)", + "ReducedJohnsonMercier(T)", + "ReducedJohnsonMercier(S)", "AlfeldSorokina(T)", "AlfeldSorokina(S)", "ArnoldQin(T, reduced=False)", diff --git a/test/FIAT/unit/test_johnson_mercier.py b/test/FIAT/unit/test_johnson_mercier.py index dcf60e6c..73d8f536 100644 --- a/test/FIAT/unit/test_johnson_mercier.py +++ b/test/FIAT/unit/test_johnson_mercier.py @@ -1,7 +1,7 @@ import pytest import numpy -from FIAT import JohnsonMercier, Nedelec +from FIAT import JohnsonMercier, ReducedJohnsonMercier, Nedelec from FIAT.reference_element import ufc_simplex from FIAT.quadrature_schemes import create_quadrature @@ -20,13 +20,16 @@ def cell(request): return K -def test_johnson_mercier_divergence_rigid_body_motions(cell): +@pytest.mark.parametrize("reduced", (False, True)) +def test_johnson_mercier_divergence_rigid_body_motions(cell, reduced): # test that the divergence of interior JM basis functions is orthogonal to # the rigid-body motions degree = 1 - variant = None sd = cell.get_spatial_dimension() - JM = JohnsonMercier(cell, degree, variant=variant) + if reduced: + JM = ReducedJohnsonMercier(cell, degree) + else: + JM = JohnsonMercier(cell, degree) ref_complex = JM.get_reference_complex() Q = create_quadrature(ref_complex, 2*(degree)-1) @@ -40,18 +43,9 @@ def test_johnson_mercier_divergence_rigid_body_motions(cell): N1_at_qpts = N1.tabulate(0, qpts) rbms = N1_at_qpts[(0,)*sd] ells = rbms * qwts[None, None, :] + L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2))) edofs = JM.entity_dofs() - idofs = edofs[sd][0] - L = numpy.tensordot(div, ells, axes=((1, 2), (1, 2))) - assert numpy.allclose(L[idofs], 0) - - if variant == "divergence": - edofs = JM.entity_dofs() - cdofs = [] - for entity in edofs[sd-1]: - cdofs.extend(edofs[sd-1][entity][:sd]) - D = L[cdofs] - M = numpy.tensordot(rbms, ells, axes=((1, 2), (1, 2))) - X = numpy.linalg.solve(M, D.T) - assert numpy.allclose(numpy.tensordot(X, rbms, axes=(0, 0)), div[cdofs]) + idofs = list(edofs[sd][0]) + L = L[idofs].reshape(len(idofs), -1) + assert numpy.allclose(L, 0) diff --git a/test/FIAT/unit/test_mtw.py b/test/FIAT/unit/test_mtw.py index c1a71120..3ca12153 100644 --- a/test/FIAT/unit/test_mtw.py +++ b/test/FIAT/unit/test_mtw.py @@ -26,7 +26,7 @@ def test_dofs(): normal_moments[bf, m] += wts[k] * nvals[bf, k] * linevals[m, k] right = np.zeros((9, 2)) right[3*ed, 0] = 1.0 - right[3*ed+2, 1] = 1.0 + right[3*ed+1, 1] = 1.0 assert np.allclose(normal_moments, right) for ed in range(3): t = T.compute_edge_tangent(ed) @@ -39,5 +39,5 @@ def test_dofs(): for k in range(len(Qline.wts)): tangent_moments[bf] += wts[k] * tvals[bf, k] * linevals[0, k] right = np.zeros(9) - right[3*ed + 1] = 1.0 + right[3*ed+2] = 1.0 assert np.allclose(tangent_moments, right) diff --git a/test/FIAT/unit/test_stokes_complex.py b/test/FIAT/unit/test_stokes_complex.py index 4e973588..a7725ed7 100644 --- a/test/FIAT/unit/test_stokes_complex.py +++ b/test/FIAT/unit/test_stokes_complex.py @@ -14,6 +14,7 @@ from FIAT.macro import CkPolynomialSet from FIAT.alfeld_sorokina import AlfeldSorokinaSpace from FIAT.arnold_qin import ArnoldQinSpace +from FIAT.bernardi_raugel import BernardiRaugel from FIAT.christiansen_hu import ChristiansenHuSpace from FIAT.guzman_neilan import ( GuzmanNeilanFirstKindH1, @@ -21,12 +22,16 @@ GuzmanNeilanH1div, GuzmanNeilanSpace) from FIAT.restricted import RestrictedElement - +from FIAT.quadrature import FacetQuadratureRule T = ufc_simplex(2) S = ufc_simplex(3) +def inner(u, v, qwts): + return numpy.tensordot(numpy.multiply(u, qwts), v, axes=(tuple(range(1, u.ndim)), )*2) + + def rHCT(cell): return RestrictedElement(HCT(cell, reduced=True), restriction_domain="vertex") @@ -182,12 +187,13 @@ def test_gn_stokes_pairs(cell, kind): check_stokes_complex(spaces, degree) +@pytest.mark.parametrize("element", (GuzmanNeilanFirstKindH1, BernardiRaugel)) @pytest.mark.parametrize("quad_scheme", (None, "KMV(2),powell-sabin")) @pytest.mark.parametrize("sd", (2, 3)) -def test_gn_dofs(sd, quad_scheme): +def test_gn_dofs(element, sd, quad_scheme): cell = symmetric_simplex(sd) - fe = GuzmanNeilanFirstKindH1(cell, 1, quad_scheme=quad_scheme) + fe = element(cell, 1, quad_scheme=quad_scheme) degree = fe.degree() assert degree == sd ref_complex = fe.get_reference_complex() @@ -206,6 +212,39 @@ def test_gn_dofs(sd, quad_scheme): assert numpy.allclose(div_moments, expected) +@pytest.mark.parametrize("sd", (2, 3)) +def test_gn_trace(sd): + degree = 1 + cell = symmetric_simplex(sd) + + gn = GuzmanNeilanFirstKindH1(cell, degree) + br = BernardiRaugel(cell, degree) + + sd = cell.get_spatial_dimension() + ref_face = cell.construct_subelement(sd-1) + Q_face = create_quadrature(ref_face, 2*degree) + + Phis = ONPolynomialSet(ref_face, degree) + phis_at_qpts = Phis.tabulate(Q_face.get_points())[(0,)*(sd-1)] + + for f in cell.topology[sd-1]: + Q = FacetQuadratureRule(cell, sd-1, f, Q_face) + vals = gn.tabulate(0, Q.get_points())[(0,)*sd] + vals -= br.tabulate(0, Q.get_points())[(0,)*sd] + + # Assert that the normal component is the same + n = cell.compute_normal(f) + normal_trace = numpy.tensordot(vals, n, axes=(1, 0)) + normal_moments = inner(normal_trace, phis_at_qpts, Q.get_weights()) + assert numpy.allclose(normal_moments, 0) + + # Assert that the tangential components are the same + for t in cell.compute_tangents(sd-1, f): + tangential_trace = numpy.tensordot(vals, t, axes=(1, 0)) + tangential_moments = inner(tangential_trace, phis_at_qpts, Q.get_weights()) + assert numpy.allclose(tangential_moments, 0) + + @pytest.mark.parametrize("cell", (T, S)) @pytest.mark.parametrize("family", ("AQ", "CH", "GN", "GN2")) def test_minimal_stokes_space(cell, family): diff --git a/test/finat/test_zany_mapping.py b/test/finat/test_zany_mapping.py index 030d2a5a..b115afde 100644 --- a/test/finat/test_zany_mapping.py +++ b/test/finat/test_zany_mapping.py @@ -150,17 +150,18 @@ def test_argyris_point(ref_to_phys): zany_piola_elements = { 2: [ - finat.MardalTaiWinther, finat.ReducedArnoldQin, finat.ArnoldWinther, finat.ArnoldWintherNC, ], 3: [ + finat.MardalTaiWinther, finat.BernardiRaugel, finat.BernardiRaugelBubble, finat.AlfeldSorokina, finat.ChristiansenHu, finat.JohnsonMercier, + finat.ReducedJohnsonMercier, finat.GuzmanNeilanFirstKindH1, finat.GuzmanNeilanSecondKindH1, finat.GuzmanNeilanBubble,