"""Types for representing a geometric domain."""
# 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
import numbers
import warnings
from ufl.cell import AbstractCell, TensorProductCell, as_cell
from ufl.core.ufl_id import attach_ufl_id
from ufl.core.ufl_type import attach_operators_from_hash_data
from ufl.corealg.traversal import traverse_unique_terminals
from ufl.finiteelement.tensorproductelement import TensorProductElement
# Export list for ufl.classes
__all_classes__ = ["AbstractDomain", "Mesh", "MeshView", "TensorProductMesh"]
[docs]class AbstractDomain(object):
"""Symbolic representation of a geometric domain with only a geometric and topological dimension."""
def __init__(self, topological_dimension, geometric_dimension):
"""Initialise."""
# Validate dimensions
if not isinstance(geometric_dimension, numbers.Integral):
raise ValueError(f"Expecting integer geometric dimension, not {geometric_dimension.__class__}")
if not isinstance(topological_dimension, numbers.Integral):
raise ValueError(f"Expecting integer topological dimension, not {topological_dimension.__class__}")
if topological_dimension > geometric_dimension:
raise ValueError("Topological dimension cannot be larger than geometric dimension.")
# Store validated dimensions
self._topological_dimension = topological_dimension
self._geometric_dimension = geometric_dimension
[docs] def geometric_dimension(self):
"""Return the dimension of the space this domain is embedded in."""
return self._geometric_dimension
[docs] def topological_dimension(self):
"""Return the dimension of the topology of this domain."""
return self._topological_dimension
# TODO: Would it be useful to have a domain representing R^d? E.g. for
# Expression.
# class EuclideanSpace(AbstractDomain):
# def __init__(self, geometric_dimension):
# AbstractDomain.__init__(self, geometric_dimension, geometric_dimension)
[docs]@attach_operators_from_hash_data
@attach_ufl_id
class Mesh(AbstractDomain):
"""Symbolic representation of a mesh."""
def __init__(self, coordinate_element, ufl_id=None, cargo=None):
"""Initialise."""
self._ufl_id = self._init_ufl_id(ufl_id)
# Store reference to object that will not be used by UFL
self._ufl_cargo = cargo
if cargo is not None and cargo.ufl_id() != self._ufl_id:
raise ValueError("Expecting cargo object (e.g. dolfin.Mesh) to have the same ufl_id.")
# No longer accepting coordinates provided as a Coefficient
from ufl.coefficient import Coefficient
if isinstance(coordinate_element, Coefficient):
raise ValueError("Expecting a coordinate element in the ufl.Mesh construct.")
# Accept a cell in place of an element for brevity Mesh(triangle)
if isinstance(coordinate_element, AbstractCell):
from ufl.finiteelement import VectorElement
cell = coordinate_element
coordinate_element = VectorElement("Lagrange", cell, 1,
dim=cell.geometric_dimension())
# Store coordinate element
self._ufl_coordinate_element = coordinate_element
# Derive dimensions from element
gdim, = coordinate_element.value_shape()
tdim = coordinate_element.cell().topological_dimension()
AbstractDomain.__init__(self, tdim, gdim)
[docs] def ufl_cargo(self):
"""Return carried object that will not be used by UFL."""
return self._ufl_cargo
[docs] def ufl_coordinate_element(self):
"""Get the coordinate element."""
return self._ufl_coordinate_element
[docs] def ufl_cell(self):
"""Get the cell."""
return self._ufl_coordinate_element.cell()
[docs] def is_piecewise_linear_simplex_domain(self):
"""Check if the domain is a piecewise linear simplex."""
return (self._ufl_coordinate_element.degree() == 1) and self.ufl_cell().is_simplex()
def __repr__(self):
"""Representation."""
r = "Mesh(%s, %s)" % (repr(self._ufl_coordinate_element), repr(self._ufl_id))
return r
def __str__(self):
"""Format as a string."""
return "<Mesh #%s>" % (self._ufl_id,)
def _ufl_hash_data_(self):
"""UFL hash data."""
return (self._ufl_id, self._ufl_coordinate_element)
def _ufl_signature_data_(self, renumbering):
"""UFL signature data."""
return ("Mesh", renumbering[self], self._ufl_coordinate_element)
# NB! Dropped __lt__ here, don't want users to write 'mesh1 <
# mesh2'.
def _ufl_sort_key_(self):
"""UFL sort key."""
typespecific = (self._ufl_id, self._ufl_coordinate_element)
return (self.geometric_dimension(), self.topological_dimension(),
"Mesh", typespecific)
[docs]@attach_operators_from_hash_data
@attach_ufl_id
class MeshView(AbstractDomain):
"""Symbolic representation of a mesh."""
def __init__(self, mesh, topological_dimension, ufl_id=None):
"""Initialise."""
self._ufl_id = self._init_ufl_id(ufl_id)
# Store mesh
self._ufl_mesh = mesh
# Derive dimensions from element
coordinate_element = mesh.ufl_coordinate_element()
gdim, = coordinate_element.value_shape()
tdim = coordinate_element.cell().topological_dimension()
AbstractDomain.__init__(self, tdim, gdim)
[docs] def ufl_mesh(self):
"""Get the mesh."""
return self._ufl_mesh
[docs] def ufl_cell(self):
"""Get the cell."""
return self._ufl_mesh.ufl_cell()
[docs] def is_piecewise_linear_simplex_domain(self):
"""Check if the domain is a piecewise linear simplex."""
return self._ufl_mesh.is_piecewise_linear_simplex_domain()
def __repr__(self):
"""Representation."""
tdim = self.topological_dimension()
r = "MeshView(%s, %s, %s)" % (repr(self._ufl_mesh), repr(tdim), repr(self._ufl_id))
return r
def __str__(self):
"""Format as a string."""
return "<MeshView #%s of dimension %d over mesh %s>" % (
self._ufl_id, self.topological_dimension(), self._ufl_mesh)
def _ufl_hash_data_(self):
"""UFL hash data."""
return (self._ufl_id,) + self._ufl_mesh._ufl_hash_data_()
def _ufl_signature_data_(self, renumbering):
"""UFL signature data."""
return ("MeshView", renumbering[self],
self._ufl_mesh._ufl_signature_data_(renumbering))
# NB! Dropped __lt__ here, don't want users to write 'mesh1 <
# mesh2'.
def _ufl_sort_key_(self):
"""UFL sort key."""
typespecific = (self._ufl_id, self._ufl_mesh)
return (self.geometric_dimension(), self.topological_dimension(),
"MeshView", typespecific)
[docs]@attach_operators_from_hash_data
@attach_ufl_id
class TensorProductMesh(AbstractDomain):
"""Symbolic representation of a mesh."""
def __init__(self, meshes, ufl_id=None):
"""Initialise."""
self._ufl_id = self._init_ufl_id(ufl_id)
# TODO: Error checking of meshes
self._ufl_meshes = meshes
# TODO: Is this what we want to do?
# Build cell from mesh cells
self._ufl_cell = TensorProductCell(*[mesh.ufl_cell() for mesh in meshes])
# TODO: Is this what we want to do?
# Build coordinate element from mesh coordinate elements
self._ufl_coordinate_element = TensorProductElement([mesh.ufl_coordinate_element() for mesh in meshes])
# Derive dimensions from meshes
gdim = sum(mesh.geometric_dimension() for mesh in meshes)
tdim = sum(mesh.topological_dimension() for mesh in meshes)
AbstractDomain.__init__(self, tdim, gdim)
[docs] def ufl_coordinate_element(self):
"""Get the coordinate element."""
return self._ufl_coordinate_element
[docs] def ufl_cell(self):
"""Get the cell."""
return self._ufl_cell
[docs] def ufl_meshes(self):
"""Get the UFL meshes."""
return self._ufl_meshes
[docs] def is_piecewise_linear_simplex_domain(self):
"""Check if the domain is a piecewise linear simplex."""
return False # TODO: Any cases this is True
def __repr__(self):
"""Representation."""
r = "TensorProductMesh(%s, %s)" % (repr(self._ufl_meshes), repr(self._ufl_id))
return r
def __str__(self):
"""Format as a string."""
return "<TensorProductMesh #%s with meshes %s>" % (
self._ufl_id, self._ufl_meshes)
def _ufl_hash_data_(self):
"""UFL hash data."""
return (self._ufl_id,) + tuple(mesh._ufl_hash_data_() for mesh in self._ufl_meshes)
def _ufl_signature_data_(self, renumbering):
"""UFL signature data."""
return ("TensorProductMesh",) + tuple(mesh._ufl_signature_data_(renumbering) for mesh in self._ufl_meshes)
# NB! Dropped __lt__ here, don't want users to write 'mesh1 <
# mesh2'.
def _ufl_sort_key_(self):
"""UFL sort key."""
typespecific = (self._ufl_id, tuple(mesh._ufl_sort_key_() for mesh in self._ufl_meshes))
return (self.geometric_dimension(), self.topological_dimension(),
"TensorProductMesh", typespecific)
# --- Utility conversion functions
[docs]def affine_mesh(cell, ufl_id=None):
"""Create a Mesh over a given cell type with an affine geometric parameterization."""
from ufl.finiteelement import VectorElement
cell = as_cell(cell)
gdim = cell.geometric_dimension()
degree = 1
coordinate_element = VectorElement("Lagrange", cell, degree, dim=gdim)
return Mesh(coordinate_element, ufl_id=ufl_id)
_default_domains = {}
[docs]def default_domain(cell):
"""Create a singular default Mesh from a cell, always returning the same Mesh object for the same cell."""
global _default_domains
warnings.warn("default_domain is deprecated.", FutureWarning)
assert isinstance(cell, AbstractCell)
domain = _default_domains.get(cell)
if domain is None:
# Create one and only one affine Mesh with a negative ufl_id
# to avoid id collision
ufl_id = -(len(_default_domains) + 1)
domain = affine_mesh(cell, ufl_id=ufl_id)
_default_domains[cell] = domain
return domain
[docs]def as_domain(domain):
"""Convert any valid object to an AbstractDomain type."""
if isinstance(domain, AbstractDomain):
# Modern UFL files and dolfin behaviour
return domain
try:
return extract_unique_domain(domain)
except AttributeError:
try:
# Legacy UFL files
# TODO: Make this conversion in the relevant constructors
# closer to the user interface?
# TODO: Make this configurable to be an error from the dolfin side?
cell = as_cell(domain)
return default_domain(cell)
except ValueError:
return domain.ufl_domain()
[docs]def sort_domains(domains):
"""Sort domains in a canonical ordering."""
return tuple(sorted(domains, key=lambda domain: domain._ufl_sort_key_()))
[docs]def join_domains(domains):
"""Take a list of domains and return a tuple with only unique domain objects.
Checks that domains with the same id are compatible.
"""
# Use hashing to join domains, ignore None
domains = set(domains) - set((None,))
if not domains:
return ()
# Check geometric dimension compatibility
gdims = set()
for domain in domains:
gdims.add(domain.geometric_dimension())
if len(gdims) != 1:
raise ValueError("Found domains with different geometric dimensions.")
gdim, = gdims
# Split into legacy and modern style domains
legacy_domains = []
modern_domains = []
for domain in domains:
if isinstance(domain, Mesh) and domain.ufl_id() < 0:
assert domain.ufl_cargo() is None
legacy_domains.append(domain)
else:
modern_domains.append(domain)
# Handle legacy domains checking
if legacy_domains:
warnings.warn("The use of Legacy domains will be deprecated by December 2023. "
"Please, use FunctionSpace instead", DeprecationWarning)
if modern_domains:
raise ValueError(
"Found both a new-style domain and a legacy default domain. "
"These should not be used interchangeably. To find the legacy "
"domain, note that it is automatically created from a cell so "
"look for constructors taking a cell.")
return tuple(legacy_domains)
# Handle modern domains checking (assuming correct by construction)
return tuple(modern_domains)
# TODO: Move these to an analysis module?
[docs]def extract_domains(expr):
"""Return all domains expression is defined on."""
domainlist = []
for t in traverse_unique_terminals(expr):
domainlist.extend(t.ufl_domains())
return sorted(join_domains(domainlist))
[docs]def extract_unique_domain(expr):
"""Return the single unique domain expression is defined on or throw an error."""
domains = extract_domains(expr)
if len(domains) == 1:
return domains[0]
elif domains:
raise ValueError("Found multiple domains, cannot return just one.")
else:
return None
[docs]def find_geometric_dimension(expr):
"""Find the geometric dimension of an expression."""
gdims = set()
for t in traverse_unique_terminals(expr):
domain = extract_unique_domain(t)
if domain is not None:
gdims.add(domain.geometric_dimension())
if hasattr(t, "ufl_element"):
element = t.ufl_element()
if element is not None:
cell = element.cell()
if cell is not None:
gdims.add(cell.geometric_dimension())
if len(gdims) != 1:
raise ValueError("Cannot determine geometric dimension from expression.")
gdim, = gdims
return gdim