From 9f1df6848b86cde1a3a8d3cda1c1c8135bffef5d Mon Sep 17 00:00:00 2001 From: Matthew Scroggs Date: Mon, 16 Oct 2023 17:24:39 +0100 Subject: [PATCH] remove ufl.legacy --- .github/workflows/tsfc-tests.yml | 2 +- ufl/legacy/__init__.py | 25 -- ufl/legacy/brokenelement.py | 53 --- ufl/legacy/elementlist.py | 481 ------------------------ ufl/legacy/enrichedelement.py | 164 --------- ufl/legacy/finiteelement.py | 235 ------------ ufl/legacy/finiteelementbase.py | 276 -------------- ufl/legacy/hdivcurl.py | 213 ----------- ufl/legacy/mixedelement.py | 562 ----------------------------- ufl/legacy/restrictedelement.py | 113 ------ ufl/legacy/tensorproductelement.py | 140 ------- 11 files changed, 1 insertion(+), 2263 deletions(-) delete mode 100644 ufl/legacy/__init__.py delete mode 100644 ufl/legacy/brokenelement.py delete mode 100644 ufl/legacy/elementlist.py delete mode 100644 ufl/legacy/enrichedelement.py delete mode 100644 ufl/legacy/finiteelement.py delete mode 100644 ufl/legacy/finiteelementbase.py delete mode 100644 ufl/legacy/hdivcurl.py delete mode 100644 ufl/legacy/mixedelement.py delete mode 100644 ufl/legacy/restrictedelement.py delete mode 100644 ufl/legacy/tensorproductelement.py diff --git a/.github/workflows/tsfc-tests.yml b/.github/workflows/tsfc-tests.yml index a861b4728..219b7a655 100644 --- a/.github/workflows/tsfc-tests.yml +++ b/.github/workflows/tsfc-tests.yml @@ -32,7 +32,7 @@ jobs: with: path: ./tsfc repository: firedrakeproject/tsfc - ref: mscroggs/newfl-legacy + ref: mscroggs/newfl-legacy2 - name: Install tsfc run: | cd tsfc diff --git a/ufl/legacy/__init__.py b/ufl/legacy/__init__.py deleted file mode 100644 index 06a0b12ca..000000000 --- a/ufl/legacy/__init__.py +++ /dev/null @@ -1,25 +0,0 @@ -"""Legacy UFL features.""" -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Andrew T. T. McRae 2014 -# Modified by Lawrence Mitchell 2014 - -import warnings as _warnings - -from ufl.legacy.brokenelement import BrokenElement -from ufl.legacy.enrichedelement import EnrichedElement, NodalEnrichedElement -from ufl.legacy.finiteelement import FiniteElement -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.legacy.hdivcurl import HCurlElement, HDivElement, WithMapping -from ufl.legacy.mixedelement import MixedElement, TensorElement, VectorElement -from ufl.legacy.restrictedelement import RestrictedElement -from ufl.legacy.tensorproductelement import TensorProductElement - -_warnings.warn("The features in ufl.legacy are deprecated and will be removed in a future version.", - FutureWarning) diff --git a/ufl/legacy/brokenelement.py b/ufl/legacy/brokenelement.py deleted file mode 100644 index 2563b868f..000000000 --- a/ufl/legacy/brokenelement.py +++ /dev/null @@ -1,53 +0,0 @@ -"""Element.""" -# -*- coding: utf-8 -*- -# Copyright (C) 2014 Andrew T. T. McRae -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Massimiliano Leoni, 2016 - -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.sobolevspace import L2 - - -class BrokenElement(FiniteElementBase): - """The discontinuous version of an existing Finite Element space.""" - def __init__(self, element): - """Init.""" - self._element = element - - family = "BrokenElement" - cell = element.cell - degree = element.degree() - quad_scheme = element.quadrature_scheme() - value_shape = element.value_shape - reference_value_shape = element.reference_value_shape - FiniteElementBase.__init__(self, family, cell, degree, - quad_scheme, value_shape, reference_value_shape) - - def __repr__(self): - """Doc.""" - return f"BrokenElement({repr(self._element)})" - - def mapping(self): - """Doc.""" - return self._element.mapping() - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - return L2 - - def reconstruct(self, **kwargs): - """Doc.""" - return BrokenElement(self._element.reconstruct(**kwargs)) - - def __str__(self): - """Doc.""" - return f"BrokenElement({repr(self._element)})" - - def shortstr(self): - """Format as string for pretty printing.""" - return f"BrokenElement({repr(self._element)})" diff --git a/ufl/legacy/elementlist.py b/ufl/legacy/elementlist.py deleted file mode 100644 index 712b94bb1..000000000 --- a/ufl/legacy/elementlist.py +++ /dev/null @@ -1,481 +0,0 @@ -"""Element. - -This module provides an extensive list of predefined finite element -families. Users or, more likely, form compilers, may register new -elements by calling the function register_element. -""" -# Copyright (C) 2008-2016 Martin Sandve Alnæs and Anders Logg -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Marie E. Rognes , 2010 -# Modified by Lizao Li , 2015, 2016 -# Modified by Massimiliano Leoni, 2016 -# Modified by Robert Kloefkorn, 2022 - -import warnings - -from numpy import asarray - -from ufl.cell import Cell, TensorProductCell -from ufl.sobolevspace import H1, H2, L2, HCurl, HDiv, HDivDiv, HEin, HInf -from ufl.utils.formatting import istr - -# List of valid elements -ufl_elements = {} - -# Aliases: aliases[name] (...) -> (standard_name, ...) -aliases = {} - - -# Function for registering new elements -def register_element(family, short_name, value_rank, sobolev_space, mapping, - degree_range, cellnames): - """Register new finite element family.""" - if family in ufl_elements: - raise ValueError(f"Finite element '{family}%s' has already been registered.") - ufl_elements[family] = (family, short_name, value_rank, sobolev_space, - mapping, degree_range, cellnames) - if short_name is not None: - ufl_elements[short_name] = (family, short_name, value_rank, sobolev_space, - mapping, degree_range, cellnames) - - -def register_alias(alias, to): - """Doc.""" - aliases[alias] = to - - -def show_elements(): - """Shows all registered elements.""" - print("Showing all registered elements:") - print("================================") - shown = set() - for k in sorted(ufl_elements.keys()): - data = ufl_elements[k] - if data in shown: - continue - shown.add(data) - (family, short_name, value_rank, sobolev_space, mapping, degree_range, cellnames) = data - print(f"Finite element family: '{family}', '{short_name}'") - print(f"Sobolev space: {sobolev_space}%s") - print(f"Mapping: {mapping}") - print(f"Degree range: {degree_range}") - print(f"Value rank: {value_rank}") - print(f"Defined on cellnames: {cellnames}") - print() - - -# FIXME: Consider cleanup of element names. Use notation from periodic -# table as the main, keep old names as compatibility aliases. - -# NOTE: Any element with polynomial degree 0 will be considered L2, -# independent of the space passed to register_element. - -# NOTE: The mapping of the element basis functions -# from reference to physical representation is -# chosen based on the sobolev space: -# HDiv = contravariant Piola, -# HCurl = covariant Piola, -# H1/L2 = no mapping. - -# TODO: If determining mapping from sobolev_space isn't sufficient in -# the future, add mapping name as another element property. - -# Cell groups -simplices = ("interval", "triangle", "tetrahedron", "pentatope") -cubes = ("interval", "quadrilateral", "hexahedron", "tesseract") -any_cell = (None, - "vertex", "interval", - "triangle", "tetrahedron", "prism", - "pyramid", "quadrilateral", "hexahedron", "pentatope", "tesseract") - -# Elements in the periodic table # TODO: Register these as aliases of -# periodic table element description instead of the other way around -register_element("Lagrange", "CG", 0, H1, "identity", (1, None), - any_cell) # "P" -register_element("Brezzi-Douglas-Marini", "BDM", 1, HDiv, - "contravariant Piola", (1, None), simplices[1:]) # "BDMF" (2d), "N2F" (3d) -register_element("Discontinuous Lagrange", "DG", 0, L2, "identity", (0, None), - any_cell) # "DP" -register_element("Discontinuous Taylor", "TDG", 0, L2, "identity", (0, None), simplices) -register_element("Nedelec 1st kind H(curl)", "N1curl", 1, HCurl, - "covariant Piola", (1, None), simplices[1:]) # "RTE" (2d), "N1E" (3d) -register_element("Nedelec 2nd kind H(curl)", "N2curl", 1, HCurl, - "covariant Piola", (1, None), simplices[1:]) # "BDME" (2d), "N2E" (3d) -register_element("Raviart-Thomas", "RT", 1, HDiv, "contravariant Piola", - (1, None), simplices[1:]) # "RTF" (2d), "N1F" (3d) - -# Elements not in the periodic table -register_element("Argyris", "ARG", 0, H2, "custom", (5, 5), ("triangle",)) -register_element("Bell", "BELL", 0, H2, "custom", (5, 5), ("triangle",)) -register_element("Brezzi-Douglas-Fortin-Marini", "BDFM", 1, HDiv, - "contravariant Piola", (1, None), simplices[1:]) -register_element("Crouzeix-Raviart", "CR", 0, L2, "identity", (1, 1), - simplices[1:]) -# TODO: Implement generic Tear operator for elements instead of this: -register_element("Discontinuous Raviart-Thomas", "DRT", 1, L2, - "contravariant Piola", (1, None), simplices[1:]) -register_element("Hermite", "HER", 0, H1, "custom", (3, 3), simplices) -register_element("Kong-Mulder-Veldhuizen", "KMV", 0, H1, "identity", (1, None), - simplices[1:]) -register_element("Mardal-Tai-Winther", "MTW", 1, H1, "contravariant Piola", (3, 3), - ("triangle",)) -register_element("Morley", "MOR", 0, H2, "custom", (2, 2), ("triangle",)) - -# Special elements -register_element("Boundary Quadrature", "BQ", 0, L2, "identity", (0, None), - any_cell) -register_element("Bubble", "B", 0, H1, "identity", (2, None), simplices) -register_element("FacetBubble", "FB", 0, H1, "identity", (2, None), simplices) -register_element("Quadrature", "Quadrature", 0, L2, "identity", (0, None), - any_cell) -register_element("Real", "R", 0, HInf, "identity", (0, 0), - any_cell + ("TensorProductCell",)) -register_element("Undefined", "U", 0, L2, "identity", (0, None), any_cell) -register_element("Radau", "Rad", 0, L2, "identity", (0, None), ("interval",)) -register_element("Regge", "Regge", 2, HEin, "double covariant Piola", - (0, None), simplices[1:]) -register_element("HDiv Trace", "HDivT", 0, L2, "identity", (0, None), any_cell) -register_element("Hellan-Herrmann-Johnson", "HHJ", 2, HDivDiv, - "double contravariant Piola", (0, None), ("triangle",)) -register_element("Nonconforming Arnold-Winther", "AWnc", 2, HDivDiv, - "double contravariant Piola", (2, 2), ("triangle", "tetrahedron")) -register_element("Conforming Arnold-Winther", "AWc", 2, HDivDiv, - "double contravariant Piola", (3, None), ("triangle", "tetrahedron")) -# Spectral elements. -register_element("Gauss-Legendre", "GL", 0, L2, "identity", (0, None), - ("interval",)) -register_element("Gauss-Lobatto-Legendre", "GLL", 0, H1, "identity", (1, None), - ("interval",)) -register_alias("Lobatto", - lambda family, dim, order, degree: ("Gauss-Lobatto-Legendre", order)) -register_alias("Lob", - lambda family, dim, order, degree: ("Gauss-Lobatto-Legendre", order)) - -register_element("Bernstein", None, 0, H1, "identity", (1, None), simplices) - - -# Let Nedelec H(div) elements be aliases to BDMs/RTs -register_alias("Nedelec 1st kind H(div)", - lambda family, dim, order, degree: ("Raviart-Thomas", order)) -register_alias("N1div", - lambda family, dim, order, degree: ("Raviart-Thomas", order)) - -register_alias("Nedelec 2nd kind H(div)", - lambda family, dim, order, degree: ("Brezzi-Douglas-Marini", - order)) -register_alias("N2div", - lambda family, dim, order, degree: ("Brezzi-Douglas-Marini", - order)) - -# Let Discontinuous Lagrange Trace element be alias to HDiv Trace -register_alias("Discontinuous Lagrange Trace", - lambda family, dim, order, degree: ("HDiv Trace", order)) -register_alias("DGT", - lambda family, dim, order, degree: ("HDiv Trace", order)) - -# New elements introduced for the periodic table 2014 -register_element("Q", None, 0, H1, "identity", (1, None), cubes) -register_element("DQ", None, 0, L2, "identity", (0, None), cubes) -register_element("RTCE", None, 1, HCurl, "covariant Piola", (1, None), - ("quadrilateral",)) -register_element("RTCF", None, 1, HDiv, "contravariant Piola", (1, None), - ("quadrilateral",)) -register_element("NCE", None, 1, HCurl, "covariant Piola", (1, None), - ("hexahedron",)) -register_element("NCF", None, 1, HDiv, "contravariant Piola", (1, None), - ("hexahedron",)) - -register_element("S", None, 0, H1, "identity", (1, None), cubes) -register_element("DPC", None, 0, L2, "identity", (0, None), cubes) -register_element("BDMCE", None, 1, HCurl, "covariant Piola", (1, None), - ("quadrilateral",)) -register_element("BDMCF", None, 1, HDiv, "contravariant Piola", (1, None), - ("quadrilateral",)) -register_element("SminusE", "SminusE", 1, HCurl, "covariant Piola", (1, None), cubes[1:3]) -register_element("SminusF", "SminusF", 1, HDiv, "contravariant Piola", (1, None), cubes[1:2]) -register_element("SminusDiv", "SminusDiv", 1, HDiv, "contravariant Piola", (1, None), cubes[1:3]) -register_element("SminusCurl", "SminusCurl", 1, HCurl, "covariant Piola", (1, None), cubes[1:3]) -register_element("AAE", None, 1, HCurl, "covariant Piola", (1, None), - ("hexahedron",)) -register_element("AAF", None, 1, HDiv, "contravariant Piola", (1, None), - ("hexahedron",)) - -# New aliases introduced for the periodic table 2014 -register_alias("P", lambda family, dim, order, degree: ("Lagrange", order)) -register_alias("DP", lambda family, dim, order, - degree: ("Discontinuous Lagrange", order)) -register_alias("RTE", lambda family, dim, order, - degree: ("Nedelec 1st kind H(curl)", order)) -register_alias("RTF", lambda family, dim, order, - degree: ("Raviart-Thomas", order)) -register_alias("N1E", lambda family, dim, order, - degree: ("Nedelec 1st kind H(curl)", order)) -register_alias("N1F", lambda family, dim, order, degree: ("Raviart-Thomas", - order)) - -register_alias("BDME", lambda family, dim, order, - degree: ("Nedelec 2nd kind H(curl)", order)) -register_alias("BDMF", lambda family, dim, order, - degree: ("Brezzi-Douglas-Marini", order)) -register_alias("N2E", lambda family, dim, order, - degree: ("Nedelec 2nd kind H(curl)", order)) -register_alias("N2F", lambda family, dim, order, - degree: ("Brezzi-Douglas-Marini", order)) - -# discontinuous elements using l2 pullbacks -register_element("DPC L2", None, 0, L2, "L2 Piola", (1, None), cubes) -register_element("DQ L2", None, 0, L2, "L2 Piola", (0, None), cubes) -register_element("Gauss-Legendre L2", "GL L2", 0, L2, "L2 Piola", (0, None), - ("interval",)) -register_element("Discontinuous Lagrange L2", "DG L2", 0, L2, "L2 Piola", (0, None), - any_cell) # "DP" - -register_alias("DP L2", lambda family, dim, order, - degree: ("Discontinuous Lagrange L2", order)) - -register_alias("P- Lambda L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) -register_alias("P Lambda L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) -register_alias("Q- Lambda L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) -register_alias("S Lambda L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) - -register_alias("P- L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) -register_alias("Q- L2", lambda family, dim, order, - degree: feec_element_l2(family, dim, order, degree)) - -# mimetic spectral elements - primal and dual complexs -register_element("Extended-Gauss-Legendre", "EGL", 0, H1, "identity", (2, None), - ("interval",)) -register_element("Extended-Gauss-Legendre Edge", "EGL-Edge", 0, L2, "identity", (1, None), - ("interval",)) -register_element("Extended-Gauss-Legendre Edge L2", "EGL-Edge L2", 0, L2, "L2 Piola", (1, None), - ("interval",)) -register_element("Gauss-Lobatto-Legendre Edge", "GLL-Edge", 0, L2, "identity", (0, None), - ("interval",)) -register_element("Gauss-Lobatto-Legendre Edge L2", "GLL-Edge L2", 0, L2, "L2 Piola", (0, None), - ("interval",)) - -# directly-defined serendipity elements ala Arbogast -# currently the theory is only really worked out for quads. -register_element("Direct Serendipity", "Sdirect", 0, H1, "physical", (1, None), - ("quadrilateral",)) -register_element("Direct Serendipity Full H(div)", "Sdirect H(div)", 1, HDiv, "physical", (1, None), - ("quadrilateral",)) -register_element("Direct Serendipity Reduced H(div)", "Sdirect H(div) red", 1, HDiv, "physical", (1, None), - ("quadrilateral",)) - - -# NOTE- the edge elements for primal mimetic spectral elements are accessed by using -# variant='mse' in the appropriate places - -def feec_element(family, n, r, k): - """Finite element exterior calculus notation. - - n = topological dimension of domain - r = polynomial order - k = form_degree - """ - # Note: We always map to edge elements in 2D, don't know how to - # differentiate otherwise? - - # Mapping from (feec name, domain dimension, form degree) to - # (family name, polynomial order) - _feec_elements = { - "P- Lambda": ( - (("P", r), ("DP", r - 1)), - (("P", r), ("RTE", r), ("DP", r - 1)), - (("P", r), ("N1E", r), ("N1F", r), ("DP", r - 1)), - ), - "P Lambda": ( - (("P", r), ("DP", r)), - (("P", r), ("BDME", r), ("DP", r)), - (("P", r), ("N2E", r), ("N2F", r), ("DP", r)), - ), - "Q- Lambda": ( - (("Q", r), ("DQ", r - 1)), - (("Q", r), ("RTCE", r), ("DQ", r - 1)), - (("Q", r), ("NCE", r), ("NCF", r), ("DQ", r - 1)), - ), - "S Lambda": ( - (("S", r), ("DPC", r)), - (("S", r), ("BDMCE", r), ("DPC", r)), - (("S", r), ("AAE", r), ("AAF", r), ("DPC", r)), - ), - } - - # New notation, old verbose notation (including "Lambda") might be - # removed - _feec_elements["P-"] = _feec_elements["P- Lambda"] - _feec_elements["P"] = _feec_elements["P Lambda"] - _feec_elements["Q-"] = _feec_elements["Q- Lambda"] - _feec_elements["S"] = _feec_elements["S Lambda"] - - family, r = _feec_elements[family][n - 1][k] - - return family, r - - -def feec_element_l2(family, n, r, k): - """Finite element exterior calculus notation. - - n = topological dimension of domain - r = polynomial order - k = form_degree - """ - # Note: We always map to edge elements in 2D, don't know how to - # differentiate otherwise? - - # Mapping from (feec name, domain dimension, form degree) to - # (family name, polynomial order) - _feec_elements = { - "P- Lambda L2": ( - (("P", r), ("DP L2", r - 1)), - (("P", r), ("RTE", r), ("DP L2", r - 1)), - (("P", r), ("N1E", r), ("N1F", r), ("DP L2", r - 1)), - ), - "P Lambda L2": ( - (("P", r), ("DP L2", r)), - (("P", r), ("BDME", r), ("DP L2", r)), - (("P", r), ("N2E", r), ("N2F", r), ("DP L2", r)), - ), - "Q- Lambda L2": ( - (("Q", r), ("DQ L2", r - 1)), - (("Q", r), ("RTCE", r), ("DQ L2", r - 1)), - (("Q", r), ("NCE", r), ("NCF", r), ("DQ L2", r - 1)), - ), - "S Lambda L2": ( - (("S", r), ("DPC L2", r)), - (("S", r), ("BDMCE", r), ("DPC L2", r)), - (("S", r), ("AAE", r), ("AAF", r), ("DPC L2", r)), - ), - } - - # New notation, old verbose notation (including "Lambda") might be - # removed - _feec_elements["P- L2"] = _feec_elements["P- Lambda L2"] - _feec_elements["P L2"] = _feec_elements["P Lambda L2"] - _feec_elements["Q- L2"] = _feec_elements["Q- Lambda L2"] - _feec_elements["S L2"] = _feec_elements["S Lambda L2"] - - family, r = _feec_elements[family][n - 1][k] - - return family, r - - -# General FEEC notation, old verbose (can be removed) -register_alias("P- Lambda", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) -register_alias("P Lambda", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) -register_alias("Q- Lambda", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) -register_alias("S Lambda", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) - -# General FEEC notation, new compact notation -register_alias("P-", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) -register_alias("Q-", lambda family, dim, order, - degree: feec_element(family, dim, order, degree)) - - -def canonical_element_description(family, cell, order, form_degree): - """Given basic element information, return corresponding element information on canonical form. - - Input: family, cell, (polynomial) order, form_degree - Output: family (canonical), short_name (for printing), order, value shape, - reference value shape, sobolev_space. - - This is used by the FiniteElement constructor to ved input - data against the element list and aliases defined in ufl. - """ - # Get domain dimensions - if cell is not None: - tdim = cell.topological_dimension() - gdim = cell.geometric_dimension() - if isinstance(cell, Cell): - cellname = cell.cellname() - else: - cellname = None - else: - tdim = None - gdim = None - cellname = None - - # Catch general FEEC notation "P" and "S" - if form_degree is not None and family in ("P", "S"): - family, order = feec_element(family, tdim, order, form_degree) - - if form_degree is not None and family in ("P L2", "S L2"): - family, order = feec_element_l2(family, tdim, order, form_degree) - - # Check whether this family is an alias for something else - while family in aliases: - if tdim is None: - raise ValueError("Need dimension to handle element aliases.") - (family, order) = aliases[family](family, tdim, order, form_degree) - - # Check that the element family exists - if family not in ufl_elements: - raise ValueError(f"Unknown finite element '{family}'.") - - # Check that element data is valid (and also get common family - # name) - (family, short_name, value_rank, sobolev_space, mapping, krange, cellnames) = ufl_elements[family] - - # Accept CG/DG on all kind of cells, but use Q/DQ on "product" cells - if cellname in set(cubes) - set(simplices) or isinstance(cell, TensorProductCell): - if family == "Lagrange": - family = "Q" - elif family == "Discontinuous Lagrange": - if order >= 1: - warnings.warn("Discontinuous Lagrange element requested on %s, creating DQ element." % cell.cellname()) - family = "DQ" - elif family == "Discontinuous Lagrange L2": - if order >= 1: - warnings.warn(f"Discontinuous Lagrange L2 element requested on {cell.cellname()}, " - "creating DQ L2 element.") - family = "DQ L2" - - # Validate cellname if a valid cell is specified - if not (cellname is None or cellname in cellnames): - raise ValueError(f"Cellname '{cellname}' invalid for '{family}' finite element.") - - # Validate order if specified - if order is not None: - if krange is None: - raise ValueError(f"Order {order} invalid for '{family}' finite element, should be None.") - kmin, kmax = krange - if not (kmin is None or (asarray(order) >= kmin).all()): - raise ValueError(f"Order {order} invalid for '{family}' finite element.") - if not (kmax is None or (asarray(order) <= kmax).all()): - raise ValueError(f"Order {istr(order)} invalid for '{family}' finite element.") - - if value_rank == 2: - # Tensor valued fundamental elements in HEin have this shape - if gdim is None or tdim is None: - raise ValueError("Cannot infer shape of element without topological and geometric dimensions.") - reference_value_shape = (tdim, tdim) - value_shape = (gdim, gdim) - elif value_rank == 1: - # Vector valued fundamental elements in HDiv and HCurl have a shape - if gdim is None or tdim is None: - raise ValueError("Cannot infer shape of element without topological and geometric dimensions.") - reference_value_shape = (tdim,) - value_shape = (gdim,) - elif value_rank == 0: - # All other elements are scalar values - reference_value_shape = () - value_shape = () - else: - raise ValueError(f"Invalid value rank {value_rank}.") - - return family, short_name, order, value_shape, reference_value_shape, sobolev_space, mapping diff --git a/ufl/legacy/enrichedelement.py b/ufl/legacy/enrichedelement.py deleted file mode 100644 index 76420f00b..000000000 --- a/ufl/legacy/enrichedelement.py +++ /dev/null @@ -1,164 +0,0 @@ -"""This module defines the UFL finite element classes.""" - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Massimiliano Leoni, 2016 - -from ufl.legacy.finiteelementbase import FiniteElementBase - - -class EnrichedElementBase(FiniteElementBase): - """The vector sum of several finite element spaces.""" - - def __init__(self, *elements): - """Doc.""" - self._elements = elements - - cell = elements[0].cell - if not all(e.cell == cell for e in elements[1:]): - raise ValueError("Cell mismatch for sub elements of enriched element.") - - if isinstance(elements[0].degree(), int): - degrees = {e.degree() for e in elements} - {None} - degree = max(degrees) if degrees else None - else: - degree = tuple(map(max, zip(*[e.degree() for e in elements]))) - - # We can allow the scheme not to be defined, but all defined - # should be equal - quad_schemes = [e.quadrature_scheme() for e in elements] - quad_schemes = [qs for qs in quad_schemes if qs is not None] - quad_scheme = quad_schemes[0] if quad_schemes else None - if not all(qs == quad_scheme for qs in quad_schemes): - raise ValueError("Quadrature scheme mismatch.") - - value_shape = elements[0].value_shape - if not all(e.value_shape == value_shape for e in elements[1:]): - raise ValueError("Element value shape mismatch.") - - reference_value_shape = elements[0].reference_value_shape - if not all(e.reference_value_shape == reference_value_shape for e in elements[1:]): - raise ValueError("Element reference value shape mismatch.") - - # mapping = elements[0].mapping() # FIXME: This fails for a mixed subelement here. - # if not all(e.mapping() == mapping for e in elements[1:]): - # raise ValueError("Element mapping mismatch.") - - # Get name of subclass: EnrichedElement or NodalEnrichedElement - class_name = self.__class__.__name__ - - # Initialize element data - FiniteElementBase.__init__(self, class_name, cell, degree, - quad_scheme, value_shape, - reference_value_shape) - - def mapping(self): - """Doc.""" - return self._elements[0].mapping() - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - elements = [e for e in self._elements] - if all(e.sobolev_space == elements[0].sobolev_space - for e in elements): - return elements[0].sobolev_space - else: - # Find smallest shared Sobolev space over all sub elements - spaces = [e.sobolev_space for e in elements] - superspaces = [{s} | set(s.parents) for s in spaces] - intersect = set.intersection(*superspaces) - for s in intersect.copy(): - for parent in s.parents: - intersect.discard(parent) - - sobolev_space, = intersect - return sobolev_space - - def variant(self): - """Doc.""" - try: - variant, = {e.variant() for e in self._elements} - return variant - except ValueError: - return None - - def reconstruct(self, **kwargs): - """Doc.""" - return type(self)(*[e.reconstruct(**kwargs) for e in self._elements]) - - @property - def embedded_subdegree(self): - """Return embedded subdegree.""" - if isinstance(self._degree, int): - return self._degree - else: - return min(e.embedded_subdegree for e in self._elements) - - @property - def embedded_superdegree(self): - """Return embedded superdegree.""" - if isinstance(self._degree, int): - return self._degree - else: - return max(e.embedded_superdegree for e in self._elements) - - -class EnrichedElement(EnrichedElementBase): - r"""The vector sum of several finite element spaces. - - .. math:: \\textrm{EnrichedElement}(V, Q) = \\{v + q | v \\in V, q \\in Q\\}. - - Dual basis is a concatenation of subelements dual bases; - primal basis is a concatenation of subelements primal bases; - resulting element is not nodal even when subelements are. - Structured basis may be exploited in form compilers. - """ - - def is_cellwise_constant(self): - """Return whether the basis functions of this element is spatially constant over each cell.""" - return all(e.is_cellwise_constant() for e in self._elements) - - def __repr__(self): - """Doc.""" - return "EnrichedElement(" + ", ".join(repr(e) for e in self._elements) + ")" - - def __str__(self): - """Format as string for pretty printing.""" - return "<%s>" % " + ".join(str(e) for e in self._elements) - - def shortstr(self): - """Format as string for pretty printing.""" - return "<%s>" % " + ".join(e.shortstr() for e in self._elements) - - -class NodalEnrichedElement(EnrichedElementBase): - r"""The vector sum of several finite element spaces. - - .. math:: \\textrm{EnrichedElement}(V, Q) = \\{v + q | v \\in V, q \\in Q\\}. - - Primal basis is reorthogonalized to dual basis which is - a concatenation of subelements dual bases; resulting - element is nodal. - """ - def is_cellwise_constant(self): - """Return whether the basis functions of this element is spatially constant over each cell.""" - return False - - def __repr__(self): - """Doc.""" - return "NodalEnrichedElement(" + ", ".join(repr(e) for e in self._elements) + ")" - - def __str__(self): - """Format as string for pretty printing.""" - return "" % ", ".join(str(e) for e in self._elements) - - def shortstr(self): - """Format as string for pretty printing.""" - return "NodalEnriched(%s)" % ", ".join(e.shortstr() for e in self._elements) diff --git a/ufl/legacy/finiteelement.py b/ufl/legacy/finiteelement.py deleted file mode 100644 index 1166ed741..000000000 --- a/ufl/legacy/finiteelement.py +++ /dev/null @@ -1,235 +0,0 @@ -"""This module defines the UFL finite element classes.""" -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Anders Logg 2014 -# Modified by Massimiliano Leoni, 2016 - -from ufl.cell import TensorProductCell, as_cell -from ufl.legacy.elementlist import canonical_element_description, simplices -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.utils.formatting import istr - - -class FiniteElement(FiniteElementBase): - """The basic finite element class for all simple finite elements.""" - # TODO: Move these to base? - __slots__ = ("_short_name", "_sobolev_space", - "_mapping", "_variant", "_repr") - - def __new__(cls, - family, - cell=None, - degree=None, - form_degree=None, - quad_scheme=None, - variant=None): - """Intercepts construction to expand CG, DG, RTCE and RTCF spaces on TensorProductCells.""" - if cell is not None: - cell = as_cell(cell) - - if isinstance(cell, TensorProductCell): - # Delay import to avoid circular dependency at module load time - from ufl.legacy.enrichedelement import EnrichedElement - from ufl.legacy.hdivcurl import HCurlElement as HCurl - from ufl.legacy.hdivcurl import HDivElement as HDiv - from ufl.legacy.tensorproductelement import TensorProductElement - - family, short_name, degree, value_shape, reference_value_shape, sobolev_space, mapping = \ - canonical_element_description(family, cell, degree, form_degree) - - if family in ["RTCF", "RTCE"]: - cell_h, cell_v = cell.sub_cells() - if cell_h.cellname() != "interval": - raise ValueError(f"{family} is available on TensorProductCell(interval, interval) only.") - if cell_v.cellname() != "interval": - raise ValueError(f"{family} is available on TensorProductCell(interval, interval) only.") - - C_elt = FiniteElement("CG", "interval", degree, variant=variant) - D_elt = FiniteElement("DG", "interval", degree - 1, variant=variant) - - CxD_elt = TensorProductElement(C_elt, D_elt, cell=cell) - DxC_elt = TensorProductElement(D_elt, C_elt, cell=cell) - - if family == "RTCF": - return EnrichedElement(HDiv(CxD_elt), HDiv(DxC_elt)) - if family == "RTCE": - return EnrichedElement(HCurl(CxD_elt), HCurl(DxC_elt)) - - elif family == "NCF": - cell_h, cell_v = cell.sub_cells() - if cell_h.cellname() != "quadrilateral": - raise ValueError(f"{family} is available on TensorProductCell(quadrilateral, interval) only.") - if cell_v.cellname() != "interval": - raise ValueError(f"{family} is available on TensorProductCell(quadrilateral, interval) only.") - - Qc_elt = FiniteElement("RTCF", "quadrilateral", degree, variant=variant) - Qd_elt = FiniteElement("DQ", "quadrilateral", degree - 1, variant=variant) - - Id_elt = FiniteElement("DG", "interval", degree - 1, variant=variant) - Ic_elt = FiniteElement("CG", "interval", degree, variant=variant) - - return EnrichedElement(HDiv(TensorProductElement(Qc_elt, Id_elt, cell=cell)), - HDiv(TensorProductElement(Qd_elt, Ic_elt, cell=cell))) - - elif family == "NCE": - cell_h, cell_v = cell.sub_cells() - if cell_h.cellname() != "quadrilateral": - raise ValueError(f"{family} is available on TensorProductCell(quadrilateral, interval) only.") - if cell_v.cellname() != "interval": - raise ValueError(f"{family} is available on TensorProductCell(quadrilateral, interval) only.") - - Qc_elt = FiniteElement("Q", "quadrilateral", degree, variant=variant) - Qd_elt = FiniteElement("RTCE", "quadrilateral", degree, variant=variant) - - Id_elt = FiniteElement("DG", "interval", degree - 1, variant=variant) - Ic_elt = FiniteElement("CG", "interval", degree, variant=variant) - - return EnrichedElement(HCurl(TensorProductElement(Qc_elt, Id_elt, cell=cell)), - HCurl(TensorProductElement(Qd_elt, Ic_elt, cell=cell))) - - elif family == "Q": - return TensorProductElement(*[FiniteElement("CG", c, degree, variant=variant) - for c in cell.sub_cells()], - cell=cell) - - elif family == "DQ": - def dq_family(cell): - """Doc.""" - return "DG" if cell.cellname() in simplices else "DQ" - return TensorProductElement(*[FiniteElement(dq_family(c), c, degree, variant=variant) - for c in cell.sub_cells()], - cell=cell) - - elif family == "DQ L2": - def dq_family_l2(cell): - """Doc.""" - return "DG L2" if cell.cellname() in simplices else "DQ L2" - return TensorProductElement(*[FiniteElement(dq_family_l2(c), c, degree, variant=variant) - for c in cell.sub_cells()], - cell=cell) - - return super(FiniteElement, cls).__new__(cls) - - def __init__(self, - family, - cell=None, - degree=None, - form_degree=None, - quad_scheme=None, - variant=None): - """Create finite element. - - Args: - family: The finite element family - cell: The geometric cell - degree: The polynomial degree (optional) - form_degree: The form degree (FEEC notation, used when field is - viewed as k-form) - quad_scheme: The quadrature scheme (optional) - variant: Hint for the local basis function variant (optional) - """ - # Note: Unfortunately, dolfin sometimes passes None for - # cell. Until this is fixed, allow it: - if cell is not None: - cell = as_cell(cell) - - ( - family, short_name, degree, value_shape, reference_value_shape, sobolev_space, mapping - ) = canonical_element_description(family, cell, degree, form_degree) - - # TODO: Move these to base? Might be better to instead - # simplify base though. - self._sobolev_space = sobolev_space - self._mapping = mapping - self._short_name = short_name or family - self._variant = variant - - # Type check variant - if variant is not None and not isinstance(variant, str): - raise ValueError("Illegal variant: must be string or None") - - # Initialize element data - FiniteElementBase.__init__(self, family, cell, degree, quad_scheme, - value_shape, reference_value_shape) - - # Cache repr string - qs = self.quadrature_scheme() - if qs is None: - quad_str = "" - else: - quad_str = ", quad_scheme=%s" % repr(qs) - v = self.variant() - if v is None: - var_str = "" - else: - var_str = ", variant=%s" % repr(v) - self._repr = "FiniteElement(%s, %s, %s%s%s)" % ( - repr(self.family()), repr(self.cell), repr(self.degree()), quad_str, var_str) - assert '"' not in self._repr - - def __repr__(self): - """Format as string for evaluation as Python object.""" - return self._repr - - def _is_globally_constant(self): - """Doc.""" - return self.family() == "Real" - - def _is_linear(self): - """Doc.""" - return self.family() == "Lagrange" and self.degree() == 1 - - def mapping(self): - """Return the mapping type for this element .""" - return self._mapping - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - return self._sobolev_space - - def variant(self): - """Return the variant used to initialise the element.""" - return self._variant - - def reconstruct(self, family=None, cell=None, degree=None, quad_scheme=None, variant=None): - """Construct a new FiniteElement object with some properties replaced with new values.""" - if family is None: - family = self.family() - if cell is None: - cell = self.cell - if degree is None: - degree = self.degree() - if quad_scheme is None: - quad_scheme = self.quadrature_scheme() - if variant is None: - variant = self.variant() - return FiniteElement(family, cell, degree, quad_scheme=quad_scheme, variant=variant) - - def __str__(self): - """Format as string for pretty printing.""" - qs = self.quadrature_scheme() - qs = "" if qs is None else "(%s)" % qs - v = self.variant() - v = "" if v is None else "(%s)" % v - return "<%s%s%s%s on a %s>" % (self._short_name, istr(self.degree()), - qs, v, self.cell) - - def shortstr(self): - """Format as string for pretty printing.""" - return f"{self._short_name}{istr(self.degree())}({self.quadrature_scheme()},{istr(self.variant())})" - - def __getnewargs__(self): - """Return the arguments which pickle needs to recreate the object.""" - return (self.family(), - self.cell, - self.degree(), - None, - self.quadrature_scheme(), - self.variant()) diff --git a/ufl/legacy/finiteelementbase.py b/ufl/legacy/finiteelementbase.py deleted file mode 100644 index 88eba6c8c..000000000 --- a/ufl/legacy/finiteelementbase.py +++ /dev/null @@ -1,276 +0,0 @@ -"""This module defines the UFL finite element classes.""" - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Massimiliano Leoni, 2016 - -from abc import abstractmethod, abstractproperty - -from ufl import pullback -from ufl.cell import AbstractCell, as_cell -from ufl.finiteelement import AbstractFiniteElement -from ufl.utils.sequences import product - - -class FiniteElementBase(AbstractFiniteElement): - """Base class for all finite elements.""" - __slots__ = ("_family", "_cell", "_degree", "_quad_scheme", - "_value_shape", "_reference_value_shape") - - # TODO: Not all these should be in the base class! In particular - # family, degree, and quad_scheme do not belong here. - def __init__(self, family, cell, degree, quad_scheme, value_shape, - reference_value_shape): - """Initialize basic finite element data.""" - if not (degree is None or isinstance(degree, (int, tuple))): - raise ValueError("Invalid degree type.") - if not isinstance(value_shape, tuple): - raise ValueError("Invalid value_shape type.") - if not isinstance(reference_value_shape, tuple): - raise ValueError("Invalid reference_value_shape type.") - - if cell is not None: - cell = as_cell(cell) - if not isinstance(cell, AbstractCell): - raise ValueError("Invalid cell type.") - - self._family = family - self._cell = cell - self._degree = degree - self._value_shape = value_shape - self._reference_value_shape = reference_value_shape - self._quad_scheme = quad_scheme - - @abstractmethod - def __repr__(self): - """Format as string for evaluation as Python object.""" - pass - - @abstractproperty - def sobolev_space(self): - """Return the underlying Sobolev space.""" - pass - - @abstractmethod - def mapping(self): - """Return the mapping type for this element.""" - pass - - def _is_globally_constant(self): - """Check if the element is a global constant. - - For Real elements, this should return True. - """ - return False - - def _is_linear(self): - """Check if the element is Lagrange degree 1.""" - return False - - def _ufl_hash_data_(self): - """Doc.""" - return repr(self) - - def _ufl_signature_data_(self): - """Doc.""" - return repr(self) - - def __hash__(self): - """Compute hash code for insertion in hashmaps.""" - return hash(self._ufl_hash_data_()) - - def __eq__(self, other): - """Compute element equality for insertion in hashmaps.""" - return type(self) is type(other) and self._ufl_hash_data_() == other._ufl_hash_data_() - - def __ne__(self, other): - """Compute element inequality for insertion in hashmaps.""" - return not self.__eq__(other) - - def __lt__(self, other): - """Compare elements by repr, to give a natural stable sorting.""" - return repr(self) < repr(other) - - def family(self): # FIXME: Undefined for base? - """Return finite element family.""" - return self._family - - def variant(self): - """Return the variant used to initialise the element.""" - return None - - def degree(self, component=None): - """Return polynomial degree of finite element.""" - return self._degree - - def quadrature_scheme(self): - """Return quadrature scheme of finite element.""" - return self._quad_scheme - - @property - def cell(self): - """Return cell of finite element.""" - return self._cell - - def is_cellwise_constant(self, component=None): - """Return whether the basis functions of this element is spatially constant over each cell.""" - return self._is_globally_constant() or self.degree() == 0 - - @property - def value_shape(self): - """Return the shape of the value space on the global domain.""" - return self._value_shape - - @property - def reference_value_shape(self): - """Return the shape of the value space on the reference cell.""" - return self._reference_value_shape - - @property - def value_size(self): - """Return the integer product of the value shape.""" - return product(self.value_shape) - - @property - def reference_value_size(self): - """Return the integer product of the reference value shape.""" - return product(self.reference_value_shape) - - def symmetry(self): # FIXME: different approach - r"""Return the symmetry dict. - - This is a mapping :math:`c_0 \\to c_1` - meaning that component :math:`c_0` is represented by component - :math:`c_1`. - A component is a tuple of one or more ints. - """ - return {} - - def _check_component(self, i): - """Check that component index i is valid.""" - sh = self.value_shape - r = len(sh) - if not (len(i) == r and all(j < k for (j, k) in zip(i, sh))): - raise ValueError( - f"Illegal component index {i} (value rank {len(i)}) " - f"for element (value rank {r}).") - - def extract_subelement_component(self, i): - """Extract direct subelement index and subelement relative component index for a given component index.""" - if isinstance(i, int): - i = (i,) - self._check_component(i) - return (None, i) - - def extract_component(self, i): - """Recursively extract component index relative to a (simple) element. - - and that element for given value component index. - """ - if isinstance(i, int): - i = (i,) - self._check_component(i) - return (i, self) - - def _check_reference_component(self, i): - """Check that reference component index i is valid.""" - sh = self.value_shape - r = len(sh) - if not (len(i) == r and all(j < k for (j, k) in zip(i, sh))): - raise ValueError( - f"Illegal component index {i} (value rank {len(i)}) " - f"for element (value rank {r}).") - - def extract_subelement_reference_component(self, i): - """Extract direct subelement index and subelement relative. - - reference component index for a given reference component index. - """ - if isinstance(i, int): - i = (i,) - self._check_reference_component(i) - return (None, i) - - def extract_reference_component(self, i): - """Recursively extract reference component index relative to a (simple) element. - - and that element for given reference value component index. - """ - if isinstance(i, int): - i = (i,) - self._check_reference_component(i) - return (i, self) - - @property - def num_sub_elements(self): - """Return number of sub-elements.""" - return 0 - - @property - def sub_elements(self): - """Return list of sub-elements.""" - return [] - - def __add__(self, other): - """Add two elements, creating an enriched element.""" - if not isinstance(other, FiniteElementBase): - raise ValueError(f"Can't add element and {other.__class__}.") - from ufl.legacy import EnrichedElement - return EnrichedElement(self, other) - - def __mul__(self, other): - """Multiply two elements, creating a mixed element.""" - if not isinstance(other, FiniteElementBase): - raise ValueError("Can't multiply element and {other.__class__}.") - from ufl.legacy import MixedElement - return MixedElement(self, other) - - def __getitem__(self, index): - """Restrict finite element to a subdomain, subcomponent or topology (cell).""" - if index in ("facet", "interior"): - from ufl.legacy import RestrictedElement - return RestrictedElement(self, index) - else: - raise KeyError(f"Invalid index for restriction: {repr(index)}") - - def __iter__(self): - """Iter.""" - raise TypeError(f"'{type(self).__name__}' object is not iterable") - - @property - def embedded_superdegree(self): - """Doc.""" - return self.degree() - - @property - def embedded_subdegree(self): - """Doc.""" - return self.degree() - - @property - def pullback(self): - """Get the pull back.""" - if self.mapping() == "identity": - return pullback.identity_pullback - elif self.mapping() == "L2 Piola": - return pullback.l2_piola - elif self.mapping() == "covariant Piola": - return pullback.covariant_piola - elif self.mapping() == "contravariant Piola": - return pullback.contravariant_piola - elif self.mapping() == "double covariant Piola": - return pullback.double_covariant_piola - elif self.mapping() == "double contravariant Piola": - return pullback.double_contravariant_piola - elif self.mapping() == "custom": - return pullback.custom_pullback - elif self.mapping() == "physical": - return pullback.physical_pullback - - raise ValueError(f"Unsupported mapping: {self.mapping()}") diff --git a/ufl/legacy/hdivcurl.py b/ufl/legacy/hdivcurl.py deleted file mode 100644 index 3f693017e..000000000 --- a/ufl/legacy/hdivcurl.py +++ /dev/null @@ -1,213 +0,0 @@ -"""Doc.""" -# Copyright (C) 2008-2016 Andrew T. T. McRae -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Massimiliano Leoni, 2016 - -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.sobolevspace import L2, HCurl, HDiv - - -class HDivElement(FiniteElementBase): - """A div-conforming version of an outer product element, assuming this makes mathematical sense.""" - __slots__ = ("_element", ) - - def __init__(self, element): - """Doc.""" - self._element = element - - family = "TensorProductElement" - cell = element.cell - degree = element.degree() - quad_scheme = element.quadrature_scheme() - value_shape = (element.cell.geometric_dimension(),) - reference_value_shape = (element.cell.topological_dimension(),) - - # Skipping TensorProductElement constructor! Bad code smell, refactor to avoid this non-inheritance somehow. - FiniteElementBase.__init__(self, family, cell, degree, - quad_scheme, value_shape, reference_value_shape) - - def __repr__(self): - """Doc.""" - return f"HDivElement({repr(self._element)})" - - def mapping(self): - """Doc.""" - return "contravariant Piola" - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - return HDiv - - def reconstruct(self, **kwargs): - """Doc.""" - return HDivElement(self._element.reconstruct(**kwargs)) - - def variant(self): - """Doc.""" - return self._element.variant() - - def __str__(self): - """Doc.""" - return f"HDivElement({repr(self._element)})" - - def shortstr(self): - """Format as string for pretty printing.""" - return f"HDivElement({self._element.shortstr()})" - - @property - def embedded_subdegree(self): - """Return embedded subdegree.""" - return self._element.embedded_subdegree - - @property - def embedded_superdegree(self): - """Return embedded superdegree.""" - return self._element.embedded_superdegree - - -class HCurlElement(FiniteElementBase): - """A curl-conforming version of an outer product element, assuming this makes mathematical sense.""" - __slots__ = ("_element",) - - def __init__(self, element): - """Doc.""" - self._element = element - - family = "TensorProductElement" - cell = element.cell - degree = element.degree() - quad_scheme = element.quadrature_scheme() - cell = element.cell - value_shape = (cell.geometric_dimension(),) - reference_value_shape = (cell.topological_dimension(),) # TODO: Is this right? - # Skipping TensorProductElement constructor! Bad code smell, - # refactor to avoid this non-inheritance somehow. - FiniteElementBase.__init__(self, family, cell, degree, quad_scheme, - value_shape, reference_value_shape) - - def __repr__(self): - """Doc.""" - return f"HCurlElement({repr(self._element)})" - - def mapping(self): - """Doc.""" - return "covariant Piola" - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - return HCurl - - def reconstruct(self, **kwargs): - """Doc.""" - return HCurlElement(self._element.reconstruct(**kwargs)) - - def variant(self): - """Doc.""" - return self._element.variant() - - def __str__(self): - """Doc.""" - return f"HCurlElement({repr(self._element)})" - - def shortstr(self): - """Format as string for pretty printing.""" - return f"HCurlElement({self._element.shortstr()})" - - -class WithMapping(FiniteElementBase): - """Specify an alternative mapping for the wrappee. - - For example, - to use identity mapping instead of Piola map with an element E, - write - remapped = WithMapping(E, "identity") - """ - - def __init__(self, wrapee, mapping): - """Doc.""" - if mapping == "symmetries": - raise ValueError("Can't change mapping to 'symmetries'") - self._mapping = mapping - self.wrapee = wrapee - - def __getattr__(self, attr): - """Doc.""" - try: - return getattr(self.wrapee, attr) - except AttributeError: - raise AttributeError("'%s' object has no attribute '%s'" % - (type(self).__name__, attr)) - - def __repr__(self): - """Doc.""" - return f"WithMapping({repr(self.wrapee)}, '{self._mapping}')" - - @property - def value_shape(self): - """Doc.""" - gdim = self.cell.geometric_dimension() - mapping = self.mapping() - if mapping in {"covariant Piola", "contravariant Piola"}: - return (gdim,) - elif mapping in {"double covariant Piola", "double contravariant Piola"}: - return (gdim, gdim) - else: - return self.wrapee.value_shape - - @property - def reference_value_shape(self): - """Doc.""" - tdim = self.cell.topological_dimension() - mapping = self.mapping() - if mapping in {"covariant Piola", "contravariant Piola"}: - return (tdim,) - elif mapping in {"double covariant Piola", "double contravariant Piola"}: - return (tdim, tdim) - else: - return self.wrapee.reference_value_shape - - def mapping(self): - """Doc.""" - return self._mapping - - @property - def sobolev_space(self): - """Return the underlying Sobolev space.""" - if self.wrapee.mapping() == self.mapping(): - return self.wrapee.sobolev_space - else: - return L2 - - def reconstruct(self, **kwargs): - """Doc.""" - mapping = kwargs.pop("mapping", self._mapping) - wrapee = self.wrapee.reconstruct(**kwargs) - return type(self)(wrapee, mapping) - - def variant(self): - """Doc.""" - return self.wrapee.variant() - - def __str__(self): - """Doc.""" - return f"WithMapping({repr(self.wrapee)}, {self._mapping})" - - def shortstr(self): - """Doc.""" - return f"WithMapping({self.wrapee.shortstr()}, {self._mapping})" - - @property - def embedded_subdegree(self): - """Return embedded subdegree.""" - return self._element.embedded_subdegree - - @property - def embedded_superdegree(self): - """Return embedded superdegree.""" - return self._element.embedded_superdegree diff --git a/ufl/legacy/mixedelement.py b/ufl/legacy/mixedelement.py deleted file mode 100644 index 97f194e55..000000000 --- a/ufl/legacy/mixedelement.py +++ /dev/null @@ -1,562 +0,0 @@ -"""This module defines the UFL finite element classes.""" -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Anders Logg 2014 -# Modified by Massimiliano Leoni, 2016 - -import numpy as np - -from ufl.cell import as_cell -from ufl.legacy.finiteelement import FiniteElement -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.permutation import compute_indices -from ufl.pullback import IdentityPullback, MixedPullback, SymmetricPullback -from ufl.utils.indexflattening import flatten_multiindex, shape_to_strides, unflatten_index -from ufl.utils.sequences import max_degree, product - - -class MixedElement(FiniteElementBase): - """A finite element composed of a nested hierarchy of mixed or simple elements.""" - __slots__ = ("_sub_elements", "_cells") - - def __init__(self, *elements, **kwargs): - """Create mixed finite element from given list of elements.""" - if type(self) is MixedElement: - if kwargs: - raise ValueError("Not expecting keyword arguments to MixedElement constructor.") - - # Un-nest arguments if we get a single argument with a list of elements - if len(elements) == 1 and isinstance(elements[0], (tuple, list)): - elements = elements[0] - # Interpret nested tuples as sub-mixedelements recursively - elements = [MixedElement(e) if isinstance(e, (tuple, list)) else e - for e in elements] - self._sub_elements = elements - - # Pick the first cell, for now all should be equal - cells = tuple(sorted(set(element.cell for element in elements) - set([None]))) - self._cells = cells - if cells: - cell = cells[0] - # Require that all elements are defined on the same cell - if not all(c == cell for c in cells[1:]): - raise ValueError("Sub elements must live on the same cell.") - else: - cell = None - - # Check that all elements use the same quadrature scheme TODO: - # We can allow the scheme not to be defined. - if len(elements) == 0: - quad_scheme = None - else: - quad_scheme = elements[0].quadrature_scheme() - if not all(e.quadrature_scheme() == quad_scheme for e in elements): - raise ValueError("Quadrature scheme mismatch for sub elements of mixed element.") - - # Compute value sizes in global and reference configurations - value_size_sum = sum(product(s.value_shape) for s in self._sub_elements) - reference_value_size_sum = sum(product(s.reference_value_shape) for s in self._sub_elements) - - # Default value shape: Treated simply as all subelement values - # unpacked in a vector. - value_shape = kwargs.get('value_shape', (value_size_sum,)) - - # Default reference value shape: Treated simply as all - # subelement reference values unpacked in a vector. - reference_value_shape = kwargs.get('reference_value_shape', (reference_value_size_sum,)) - - # Validate value_shape (deliberately not for subclasses - # VectorElement and TensorElement) - if type(self) is MixedElement: - # This is not valid for tensor elements with symmetries, - # assume subclasses deal with their own validation - if product(value_shape) != value_size_sum: - raise ValueError("Provided value_shape doesn't match the " - "total value size of all subelements.") - - # Initialize element data - degrees = {e.degree() for e in self._sub_elements} - {None} - degree = max_degree(degrees) if degrees else None - FiniteElementBase.__init__(self, "Mixed", cell, degree, quad_scheme, - value_shape, reference_value_shape) - - def __repr__(self): - """Doc.""" - return "MixedElement(" + ", ".join(repr(e) for e in self._sub_elements) + ")" - - def _is_linear(self): - """Doc.""" - return all(i._is_linear() for i in self._sub_elements) - - def reconstruct_from_elements(self, *elements): - """Reconstruct a mixed element from new subelements.""" - if all(a == b for (a, b) in zip(elements, self._sub_elements)): - return self - return MixedElement(*elements) - - def symmetry(self): - r"""Return the symmetry dict, which is a mapping :math:`c_0 \\to c_1`. - - meaning that component :math:`c_0` is represented by component - :math:`c_1`. - A component is a tuple of one or more ints. - """ - # Build symmetry map from symmetries of subelements - sm = {} - # Base index of the current subelement into mixed value - j = 0 - for e in self._sub_elements: - sh = e.value_shape - st = shape_to_strides(sh) - # Map symmetries of subelement into index space of this - # element - for c0, c1 in e.symmetry().items(): - j0 = flatten_multiindex(c0, st) + j - j1 = flatten_multiindex(c1, st) + j - sm[(j0,)] = (j1,) - # Update base index for next element - j += product(sh) - if j != product(self.value_shape): - raise ValueError("Size mismatch in symmetry algorithm.") - return sm or {} - - @property - def sobolev_space(self): - """Doc.""" - return max(e.sobolev_space for e in self._sub_elements) - - def mapping(self): - """Doc.""" - if all(e.mapping() == "identity" for e in self._sub_elements): - return "identity" - else: - return "undefined" - - @property - def num_sub_elements(self): - """Return number of sub elements.""" - return len(self._sub_elements) - - @property - def sub_elements(self): - """Return list of sub elements.""" - return self._sub_elements - - def extract_subelement_component(self, i): - """Extract direct subelement index and subelement relative. - - component index for a given component index. - """ - if isinstance(i, int): - i = (i,) - self._check_component(i) - - # Select between indexing modes - if len(self.value_shape) == 1: - # Indexing into a long vector of flattened subelement - # shapes - j, = i - - # Find subelement for this index - for sub_element_index, e in enumerate(self._sub_elements): - sh = e.value_shape - si = product(sh) - if j < si: - break - j -= si - if j < 0: - raise ValueError("Moved past last value component!") - - # Convert index into a shape tuple - st = shape_to_strides(sh) - component = unflatten_index(j, st) - else: - # Indexing into a multidimensional tensor where subelement - # index is first axis - sub_element_index = i[0] - if sub_element_index >= len(self._sub_elements): - raise ValueError(f"Illegal component index (dimension {sub_element_index}).") - component = i[1:] - return (sub_element_index, component) - - def extract_component(self, i): - """Recursively extract component index relative to a (simple) element. - - and that element for given value component index. - """ - sub_element_index, component = self.extract_subelement_component(i) - return self._sub_elements[sub_element_index].extract_component(component) - - def extract_subelement_reference_component(self, i): - """Extract direct subelement index and subelement relative. - - reference_component index for a given reference_component index. - """ - if isinstance(i, int): - i = (i,) - self._check_reference_component(i) - - # Select between indexing modes - assert len(self.reference_value_shape) == 1 - # Indexing into a long vector of flattened subelement shapes - j, = i - - # Find subelement for this index - for sub_element_index, e in enumerate(self._sub_elements): - sh = e.reference_value_shape - si = product(sh) - if j < si: - break - j -= si - if j < 0: - raise ValueError("Moved past last value reference_component!") - - # Convert index into a shape tuple - st = shape_to_strides(sh) - reference_component = unflatten_index(j, st) - return (sub_element_index, reference_component) - - def extract_reference_component(self, i): - """Recursively extract reference_component index relative to a (simple) element. - - and that element for given value reference_component index. - """ - sub_element_index, reference_component = self.extract_subelement_reference_component(i) - return self._sub_elements[sub_element_index].extract_reference_component(reference_component) - - def is_cellwise_constant(self, component=None): - """Return whether the basis functions of this element is spatially constant over each cell.""" - if component is None: - return all(e.is_cellwise_constant() for e in self.sub_elements) - else: - i, e = self.extract_component(component) - return e.is_cellwise_constant() - - def degree(self, component=None): - """Return polynomial degree of finite element.""" - if component is None: - return self._degree # from FiniteElementBase, computed as max of subelements in __init__ - else: - i, e = self.extract_component(component) - return e.degree() - - @property - def embedded_subdegree(self): - """Return embedded subdegree.""" - if isinstance(self._degree, int): - return self._degree - else: - return min(e.embedded_subdegree for e in self.sub_elements) - - @property - def embedded_superdegree(self): - """Return embedded superdegree.""" - if isinstance(self._degree, int): - return self._degree - else: - return max(e.embedded_superdegree for e in self.sub_elements) - - def reconstruct(self, **kwargs): - """Doc.""" - return MixedElement(*[e.reconstruct(**kwargs) for e in self.sub_elements]) - - def variant(self): - """Doc.""" - try: - variant, = {e.variant() for e in self.sub_elements} - return variant - except ValueError: - return None - - def __str__(self): - """Format as string for pretty printing.""" - tmp = ", ".join(str(element) for element in self._sub_elements) - return "" - - def shortstr(self): - """Format as string for pretty printing.""" - tmp = ", ".join(element.shortstr() for element in self._sub_elements) - return "Mixed<" + tmp + ">" - - @property - def pullback(self): - """Get the pull back.""" - for e in self.sub_elements: - if not isinstance(e.pullback, IdentityPullback): - return MixedPullback(self) - return IdentityPullback() - - -class VectorElement(MixedElement): - """A special case of a mixed finite element where all elements are equal.""" - - __slots__ = ("_repr", "_mapping", "_sub_element") - - def __init__(self, family, cell=None, degree=None, dim=None, - form_degree=None, quad_scheme=None, variant=None): - """Create vector element (repeated mixed element).""" - if isinstance(family, FiniteElementBase): - sub_element = family - cell = sub_element.cell - variant = sub_element.variant() - else: - if cell is not None: - cell = as_cell(cell) - # Create sub element - sub_element = FiniteElement(family, cell, degree, - form_degree=form_degree, - quad_scheme=quad_scheme, - variant=variant) - - # Set default size if not specified - if dim is None: - if cell is None: - raise ValueError("Cannot infer vector dimension without a cell.") - dim = cell.geometric_dimension() - - self._mapping = sub_element.mapping() - # Create list of sub elements for mixed element constructor - sub_elements = [sub_element] * dim - - # Compute value shapes - value_shape = (dim,) + sub_element.value_shape - reference_value_shape = (dim,) + sub_element.reference_value_shape - - # Initialize element data - MixedElement.__init__(self, sub_elements, value_shape=value_shape, - reference_value_shape=reference_value_shape) - - FiniteElementBase.__init__(self, sub_element.family(), sub_element.cell, sub_element.degree(), - sub_element.quadrature_scheme(), value_shape, reference_value_shape) - - self._sub_element = sub_element - - if variant is None: - var_str = "" - else: - var_str = ", variant='" + variant + "'" - - # Cache repr string - self._repr = f"VectorElement({repr(sub_element)}, dim={dim}{var_str})" - - def __repr__(self): - """Doc.""" - return self._repr - - def reconstruct(self, **kwargs): - """Doc.""" - sub_element = self._sub_element.reconstruct(**kwargs) - return VectorElement(sub_element, dim=len(self.sub_elements)) - - def variant(self): - """Return the variant used to initialise the element.""" - return self._sub_element.variant() - - def mapping(self): - """Doc.""" - return self._mapping - - def __str__(self): - """Format as string for pretty printing.""" - return ("" % - (len(self._sub_elements), self._sub_element)) - - def shortstr(self): - """Format as string for pretty printing.""" - return "Vector<%d x %s>" % (len(self._sub_elements), - self._sub_element.shortstr()) - - -class TensorElement(MixedElement): - """A special case of a mixed finite element where all elements are equal.""" - __slots__ = ("_sub_element", "_shape", "_symmetry", - "_sub_element_mapping", - "_flattened_sub_element_mapping", - "_mapping", "_repr") - - def __init__(self, family, cell=None, degree=None, shape=None, - symmetry=None, quad_scheme=None, variant=None): - """Create tensor element (repeated mixed element with optional symmetries).""" - if isinstance(family, FiniteElementBase): - sub_element = family - cell = sub_element.cell - variant = sub_element.variant() - else: - if cell is not None: - cell = as_cell(cell) - # Create scalar sub element - sub_element = FiniteElement(family, cell, degree, quad_scheme=quad_scheme, - variant=variant) - - # Set default shape if not specified - if shape is None: - if cell is None: - raise ValueError("Cannot infer tensor shape without a cell.") - dim = cell.geometric_dimension() - shape = (dim, dim) - - if symmetry is None: - symmetry = {} - elif symmetry is True: - # Construct default symmetry dict for matrix elements - if not (len(shape) == 2 and shape[0] == shape[1]): - raise ValueError("Cannot set automatic symmetry for non-square tensor.") - symmetry = dict(((i, j), (j, i)) for i in range(shape[0]) - for j in range(shape[1]) if i > j) - else: - if not isinstance(symmetry, dict): - raise ValueError("Expecting symmetry to be None (unset), True, or dict.") - - # Validate indices in symmetry dict - for i, j in symmetry.items(): - if len(i) != len(j): - raise ValueError("Non-matching length of symmetry index tuples.") - for k in range(len(i)): - if not (i[k] >= 0 and j[k] >= 0 and i[k] < shape[k] and j[k] < shape[k]): - raise ValueError("Symmetry dimensions out of bounds.") - - # Compute all index combinations for given shape - indices = compute_indices(shape) - - # Compute mapping from indices to sub element number, - # accounting for symmetry - sub_elements = [] - sub_element_mapping = {} - for index in indices: - if index in symmetry: - continue - sub_element_mapping[index] = len(sub_elements) - sub_elements += [sub_element] - - # Update mapping for symmetry - for index in indices: - if index in symmetry: - sub_element_mapping[index] = sub_element_mapping[symmetry[index]] - flattened_sub_element_mapping = [sub_element_mapping[index] for i, - index in enumerate(indices)] - - # Compute value shape - value_shape = shape - - # Compute reference value shape based on symmetries - if symmetry: - reference_value_shape = (product(shape) - len(symmetry),) - self._mapping = "symmetries" - else: - reference_value_shape = shape - self._mapping = sub_element.mapping() - - value_shape = value_shape + sub_element.value_shape - reference_value_shape = reference_value_shape + sub_element.reference_value_shape - # Initialize element data - MixedElement.__init__(self, sub_elements, value_shape=value_shape, - reference_value_shape=reference_value_shape) - self._family = sub_element.family() - self._degree = sub_element.degree() - self._sub_element = sub_element - self._shape = shape - self._symmetry = symmetry - self._sub_element_mapping = sub_element_mapping - self._flattened_sub_element_mapping = flattened_sub_element_mapping - - if variant is None: - var_str = "" - else: - var_str = ", variant='" + variant + "'" - - # Cache repr string - self._repr = (f"TensorElement({repr(sub_element)}, shape={shape}, " - f"symmetry={symmetry}{var_str})") - - @property - def pullback(self): - """Get pull back.""" - if len(self._symmetry) > 0: - sub_element_value_shape = self.sub_elements[0].value_shape - for e in self.sub_elements: - if e.value_shape != sub_element_value_shape: - raise ValueError("Sub-elements must all have the same value size") - symmetry = {} - n = 0 - for i in np.ndindex(self.value_shape[:len(self.value_shape)-len(sub_element_value_shape)]): - if i in self._symmetry and self._symmetry[i] in symmetry: - symmetry[i] = symmetry[self._symmetry[i]] - else: - symmetry[i] = n - n += 1 - return SymmetricPullback(self, symmetry) - return super().pullback - - def __repr__(self): - """Doc.""" - return self._repr - - def variant(self): - """Return the variant used to initialise the element.""" - return self._sub_element.variant() - - def mapping(self): - """Doc.""" - return self._mapping - - def flattened_sub_element_mapping(self): - """Doc.""" - return self._flattened_sub_element_mapping - - def extract_subelement_component(self, i): - """Extract direct subelement index and subelement relative. - - component index for a given component index. - """ - if isinstance(i, int): - i = (i,) - self._check_component(i) - - i = self.symmetry().get(i, i) - l = len(self._shape) # noqa: E741 - ii = i[:l] - jj = i[l:] - if ii not in self._sub_element_mapping: - raise ValueError(f"Illegal component index {i}.") - k = self._sub_element_mapping[ii] - return (k, jj) - - def symmetry(self): - r"""Return the symmetry dict, which is a mapping :math:`c_0 \\to c_1`. - - meaning that component :math:`c_0` is represented by component - :math:`c_1`. - A component is a tuple of one or more ints. - """ - return self._symmetry - - def reconstruct(self, **kwargs): - """Doc.""" - sub_element = self._sub_element.reconstruct(**kwargs) - return TensorElement(sub_element, shape=self._shape, symmetry=self._symmetry) - - def __str__(self): - """Format as string for pretty printing.""" - if self._symmetry: - tmp = ", ".join("%s -> %s" % (a, b) for (a, b) in self._symmetry.items()) - sym = " with symmetries (%s)" % tmp - else: - sym = "" - return ("" % - (self.value_shape, self._sub_element, sym)) - - def shortstr(self): - """Format as string for pretty printing.""" - if self._symmetry: - tmp = ", ".join("%s -> %s" % (a, b) for (a, b) in self._symmetry.items()) - sym = " with symmetries (%s)" % tmp - else: - sym = "" - return "Tensor<%s x %s%s>" % (self.value_shape, - self._sub_element.shortstr(), sym) diff --git a/ufl/legacy/restrictedelement.py b/ufl/legacy/restrictedelement.py deleted file mode 100644 index 767444087..000000000 --- a/ufl/legacy/restrictedelement.py +++ /dev/null @@ -1,113 +0,0 @@ -"""This module defines the UFL finite element classes.""" - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Massimiliano Leoni, 2016 - -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.sobolevspace import L2 - -valid_restriction_domains = ("interior", "facet", "face", "edge", "vertex") - - -class RestrictedElement(FiniteElementBase): - """Represents the restriction of a finite element to a type of cell entity.""" - - def __init__(self, element, restriction_domain): - """Doc.""" - if not isinstance(element, FiniteElementBase): - raise ValueError("Expecting a finite element instance.") - if restriction_domain not in valid_restriction_domains: - raise ValueError(f"Expecting one of the strings: {valid_restriction_domains}") - - FiniteElementBase.__init__(self, "RestrictedElement", element.cell, - element.degree(), - element.quadrature_scheme(), - element.value_shape, - element.reference_value_shape) - - self._element = element - - self._restriction_domain = restriction_domain - - def __repr__(self): - """Doc.""" - return f"RestrictedElement({repr(self._element)}, {repr(self._restriction_domain)})" - - @property - def sobolev_space(self): - """Doc.""" - if self._restriction_domain == "interior": - return L2 - else: - return self._element.sobolev_space - - def is_cellwise_constant(self): - """Return whether the basis functions of this element is spatially constant over each cell.""" - return self._element.is_cellwise_constant() - - def _is_linear(self): - """Doc.""" - return self._element._is_linear() - - def sub_element(self): - """Return the element which is restricted.""" - return self._element - - def mapping(self): - """Doc.""" - return self._element.mapping() - - def restriction_domain(self): - """Return the domain onto which the element is restricted.""" - return self._restriction_domain - - def reconstruct(self, **kwargs): - """Doc.""" - element = self._element.reconstruct(**kwargs) - return RestrictedElement(element, self._restriction_domain) - - def __str__(self): - """Format as string for pretty printing.""" - return "<%s>|_{%s}" % (self._element, self._restriction_domain) - - def shortstr(self): - """Format as string for pretty printing.""" - return "<%s>|_{%s}" % (self._element.shortstr(), - self._restriction_domain) - - def symmetry(self): - r"""Return the symmetry dict, which is a mapping :math:`c_0 \\to c_1`. - - meaning that component :math:`c_0` is represented by component - :math:`c_1`. A component is a tuple of one or more ints. - """ - return self._element.symmetry() - - @property - def num_sub_elements(self): - """Return number of sub elements.""" - return self._element.num_sub_elements - - @property - def sub_elements(self): - """Return list of sub elements.""" - return self._element.sub_elements - - def num_restricted_sub_elements(self): - """Return number of restricted sub elements.""" - return 1 - - def restricted_sub_elements(self): - """Return list of restricted sub elements.""" - return (self._element,) - - def variant(self): - """Doc.""" - return self._element.variant() diff --git a/ufl/legacy/tensorproductelement.py b/ufl/legacy/tensorproductelement.py deleted file mode 100644 index d9f1ff1bd..000000000 --- a/ufl/legacy/tensorproductelement.py +++ /dev/null @@ -1,140 +0,0 @@ -"""This module defines the UFL finite element classes.""" - -# Copyright (C) 2008-2016 Martin Sandve Alnæs -# -# This file is part of UFL (https://www.fenicsproject.org) -# -# SPDX-License-Identifier: LGPL-3.0-or-later -# -# Modified by Kristian B. Oelgaard -# Modified by Marie E. Rognes 2010, 2012 -# Modified by Massimiliano Leoni, 2016 - -from itertools import chain - -from ufl.cell import TensorProductCell, as_cell -from ufl.legacy.finiteelementbase import FiniteElementBase -from ufl.sobolevspace import DirectionalSobolevSpace - - -class TensorProductElement(FiniteElementBase): - r"""The tensor product of :math:`d` element spaces. - - .. math:: V = V_1 \otimes V_2 \otimes ... \otimes V_d - - Given bases :math:`\{\phi_{j_i}\}` of the spaces :math:`V_i` for :math:`i = 1, ...., d`, - :math:`\{ \phi_{j_1} \otimes \phi_{j_2} \otimes \cdots \otimes \phi_{j_d} - \}` forms a basis for :math:`V`. - """ - __slots__ = ("_sub_elements", "_cell") - - def __init__(self, *elements, **kwargs): - """Create TensorProductElement from a given list of elements.""" - if not elements: - raise ValueError("Cannot create TensorProductElement from empty list.") - - keywords = list(kwargs.keys()) - if keywords and keywords != ["cell"]: - raise ValueError("TensorProductElement got an unexpected keyword argument '%s'" % keywords[0]) - cell = kwargs.get("cell") - - family = "TensorProductElement" - - if cell is None: - # Define cell as the product of each elements cell - cell = TensorProductCell(*[e.cell for e in elements]) - else: - cell = as_cell(cell) - - # Define polynomial degree as a tuple of sub-degrees - degree = tuple(e.degree() for e in elements) - - # No quadrature scheme defined - quad_scheme = None - - # match FIAT implementation - value_shape = tuple(chain(*[e.value_shape for e in elements])) - reference_value_shape = tuple(chain(*[e.reference_value_shape for e in elements])) - if len(value_shape) > 1: - raise ValueError("Product of vector-valued elements not supported") - if len(reference_value_shape) > 1: - raise ValueError("Product of vector-valued elements not supported") - - FiniteElementBase.__init__(self, family, cell, degree, - quad_scheme, value_shape, - reference_value_shape) - self._sub_elements = elements - self._cell = cell - - def __repr__(self): - """Doc.""" - return "TensorProductElement(" + ", ".join(repr(e) for e in self._sub_elements) + f", cell={repr(self._cell)})" - - def mapping(self): - """Doc.""" - if all(e.mapping() == "identity" for e in self._sub_elements): - return "identity" - elif all(e.mapping() == "L2 Piola" for e in self._sub_elements): - return "L2 Piola" - else: - return "undefined" - - @property - def sobolev_space(self): - """Return the underlying Sobolev space of the TensorProductElement.""" - elements = self._sub_elements - if all(e.sobolev_space == elements[0].sobolev_space - for e in elements): - return elements[0].sobolev_space - else: - # Generate a DirectionalSobolevSpace which contains - # continuity information parametrized by spatial index - orders = [] - for e in elements: - e_dim = e.cell.geometric_dimension() - e_order = (e.sobolev_space._order,) * e_dim - orders.extend(e_order) - return DirectionalSobolevSpace(orders) - - @property - def num_sub_elements(self): - """Return number of subelements.""" - return len(self._sub_elements) - - @property - def sub_elements(self): - """Return subelements (factors).""" - return self._sub_elements - - def reconstruct(self, **kwargs): - """Doc.""" - cell = kwargs.pop("cell", self.cell) - return TensorProductElement(*[e.reconstruct(**kwargs) for e in self.sub_elements], cell=cell) - - def variant(self): - """Doc.""" - try: - variant, = {e.variant() for e in self.sub_elements} - return variant - except ValueError: - return None - - def __str__(self): - """Pretty-print.""" - return "TensorProductElement(%s, cell=%s)" \ - % (', '.join([str(e) for e in self._sub_elements]), str(self._cell)) - - def shortstr(self): - """Short pretty-print.""" - return "TensorProductElement(%s, cell=%s)" \ - % (', '.join([e.shortstr() for e in self._sub_elements]), str(self._cell)) - - @property - def embedded_superdegree(self): - """Doc.""" - return sum(self.degree()) - - @property - def embedded_subdegree(self): - """Doc.""" - return min(self.degree())