"""
Cell domain descriptions
"""
import abc
from collections import namedtuple
from typing import Iterable, NamedTuple, Optional
import dolfin as fem
import numpy as np
from ufl.measure import Measure # type: ignore
from mtnlion.domain import Domain
from mtnlion.variable import UndefinedDomain
[docs]class DomainData(NamedTuple):
"""
Data required to minimally describe the problem domain.
"""
mesh: fem.cpp.mesh.Mesh # pylint: disable=no-member
domain_measure: fem.Measure
boundary_measure: fem.Measure
[docs]class MeasureMap(NamedTuple):
"""
A mapping between measures, multiplicities, and primary (parent) domains
"""
measure: Optional[fem.Measure]
multiple: int
primary_domain: Optional[str]
[docs]class DomainInterface(abc.ABC):
"""
Minimal interface for describing problem domains.
"""
BoundaryMap = namedtuple("BoundaryMap", "primary_domain, marker")
def __init__(self, mesh_grid: np.ndarray):
"""
Create a `DomainInterface` object.
:param mesh_grid: numpy array describing the shape of the mesh
"""
self.mesh_grid = mesh_grid
self.domains: Domain[str, MeasureMap] = Domain({})
self.domain_data = self._domain_data()
@abc.abstractmethod
def _domain_data(self) -> DomainData:
"""
Override this method to create a `DomainData` object.
:return: `DomainData`
"""
[docs] def add_domain_measure(self, domain: str, measure: Measure, multiple: int = 1) -> None:
"""
Set the domain measure and its multiplicity for a given domain.
:param domain: The domain to assign
:param measure: The measure to apply to the domain
:param multiple: The multiplicity to apply to the domain measure
:return: None
"""
update = self.domains.get(domain, [])
update += [MeasureMap(measure, multiple, None)]
self.domains[domain] = update
[docs] def boundary_of(
self,
boundary_domain: str,
primary_domains: Iterable[str],
measures: Iterable[Measure],
multiples: Iterable[int],
) -> None:
"""
Create a boundary domain that is the boundary of one or more primary domains and assign a boundary measure
with its corresponding multiplicity.
:param boundary_domain: name of the boundary domain
:param primary_domains: names of the domains this domain is a boundary of
:param measures: assign a boundary measure to the boundary domain
:param multiples: assign a multiplicity to the boundary measure
:return: None
"""
update = self.domains.get(boundary_domain, [])
update += [MeasureMap(dm, ml, pd) for pd, ml, dm in zip(primary_domains, multiples, measures)]
self.domains[boundary_domain] = update
[docs] def is_boundary(self, domain: str) -> bool:
"""
Deduce if a given domain is a boundary domain by whether or not is has parent domains defined.
:param domain: name of the domain to check
:return: true if domain is a boundary
"""
if domain in self.domains:
return any(map(lambda x: isinstance(x.primary_domain, str), self.domains[domain]))
raise KeyError(repr(UndefinedDomain(domain)))
[docs] def create_element(self, element_type: str, degree: int) -> fem.FiniteElement:
"""
Create a FEniCS finite element. See documentation for `fem.FiniteElement`.
:param element_type: Type of the finite element.
:param degree: Degree of the element
:return: `fem.FiniteElement`
"""
return fem.FiniteElement(element_type, self.domain_data.mesh.ufl_cell(), degree)
[docs] def create_function_space(self, element: fem.FiniteElement):
"""
Create a FEniCS function space given an element. See documentation for `fem.FunctionSpace`.
:param element: Element to use to create a function space
:return: `fem.FunctionSpace`
"""
return fem.FunctionSpace(self.domain_data.mesh, element)
[docs]class P2D(DomainInterface):
"""
Definition for a pseudo two-dimensional domain.
"""
def __init__(self, mesh_grid: np.ndarray):
"""
Create a P2D domain with nodes specified at the given points.
:param mesh_grid: Points at which the one-dimensional domain should be evaluated
"""
super(P2D, self).__init__(mesh_grid)
for domain in ["anode", "separator", "cathode"]:
self.add_domain_measure(domain, self.domain_data.domain_measure)
# anode_cc is the left boundary of anode
self.boundary_of("anode_cc", ["anode"], [self.domain_data.boundary_measure(0)], [1])
# anode-separator is the right boundary of anode and the left boundary of separator
self.boundary_of(
"anode-separator",
["anode", "separator"],
[self.domain_data.boundary_measure(1), self.domain_data.boundary_measure(0)],
[1, -1],
)
# separator-cathode is the right boundary of separator and the left boundary of cathode
self.boundary_of(
"separator-cathode",
["separator", "cathode"],
[self.domain_data.boundary_measure(1), self.domain_data.boundary_measure(0)],
[-1, 1],
)
# cathode_cc is the right boundary of cathode
self.boundary_of("cathode_cc", ["cathode"], [self.domain_data.boundary_measure(1)], [1])
def _domain_data(self) -> DomainData:
"""
Create a `DomainData` object.
:return: `DomainData`
"""
mesh = fem.IntervalMesh(len(self.mesh_grid) - 1, 0, 3)
mesh.coordinates()[:] = np.array([self.mesh_grid]).transpose()
class Left(fem.SubDomain):
"""
Define the left boundary of the domain.
"""
# pylint: disable=too-few-public-methods
# pylint: disable=invalid-name
def inside(self, x, on_boundary):
# pylint: disable=missing-function-docstring
# pylint: disable=no-self-use
return x[0] < 0.0 + fem.DOLFIN_EPS and on_boundary
class Right(fem.SubDomain):
"""
Define the right boundary of the domain.
"""
# pylint: disable=too-few-public-methods
# pylint: disable=invalid-name
def inside(self, x, on_boundary):
# pylint: disable=missing-function-docstring
# pylint: disable=no-self-use
return x[0] > 1.0 - fem.DOLFIN_EPS and on_boundary
boundaries = fem.MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(2) # pylint: disable=no-member
left = Left()
right = Right()
left.mark(boundaries, 0)
right.mark(boundaries, 1)
domain_measure = fem.Measure("dx", domain=mesh)
boundary_measure = fem.Measure("ds", domain=mesh, subdomain_data=boundaries)
return DomainData(mesh, domain_measure, boundary_measure)