from collections import defaultdict
import logging
import warnings
import numpy as np
from cached_property import cached_property
from sortedcontainers import SortedList
import dolfin
from ..util import SetattrInitMixin
from .delayed_form import DelayedForm
from .expr import constantify_if_literal
logger = logging.getLogger(__name__)
[docs]class BoundaryCondition(SetattrInitMixin):
'''Boundary condition class, used both for essential and natural BCs.
Parameters
----------
priority:
Priority of this boundary condition. Lower values have higher
priority.
region: :py:class:`.topology.FacetRegion`
Facet region where this BC applies.
value: :py:class:`ufl.core.expr.Expr`
Value of the boundary condition.
require_natural: bool, optional
Does this BC require natural (as opposed to inflexible essential)
application? By default false.
internal: bool
Apply this BC to internal facets? (default: False)
external: bool
Apply this BC to external facets? (default: True)
should_extract: bool
Whether to automatically apply extractor to this BC.
'''
require_natural = False
internal = False
external = True
should_extract = True
[docs]class ValueRule(SetattrInitMixin):
'''Rule for spatially-dependent value.
Parameters
----------
priority:
Priority of this rule. Lower values have higher priority.
region: :py:class:`.topology.CellRegion`
Cell region where this value applies.
value: :py:class:`ufl.core.expr.Expr`
Value.
'''
[docs]class CombinedBoundaryCondition(SetattrInitMixin):
'''Combination of multiple :py:class:`BoundaryCondition`.
Parameters
----------
bcs: SortedList
Boundary conditions to combine.
value_extractor: callable, optional
Applied to :py:class:`BoundaryCondition` objects to obtain the
boundary condition value. See :py:attr:`map`.
mesh_util: :py:class:`.mesh_util.MeshUtil`
Used to evaluate abstract topology into actual facet values, and
other things.
essential_var: :py:class:`dolfin.Function`, optional
Variable to apply BC on (used for essential BCs).
natural_ignore_sign: bool, optional
Whether to ignore the facet sign when building a natural BC.
Attributes
----------
map: dict
Mapping :code:`{(facet_value, sign): extracted_bc}`.
reverse_map: dict
Mapping :code:`{extracted_bc: {(facet_value, sign), ...}}`.
fv_set: set
Set of referenced facet values. Useful to check that BCs don't
overlap.
'''
value_extractor = None
@property
def mesh_data(self):
return self.mesh_util.mesh_data
@cached_property
def map(self):
bcs = self.bcs
func = self.value_extractor
if func is None:
func = lambda x: x.value
mapping = {}
for bc in bcs:
region = bc.region
fvss = self.mesh_data.evaluate_topology(region)
for fvs in fvss:
flipped = (fvs[0], -fvs[1])
if flipped not in mapping:
extracted_bc = ExtractedBoundaryCondition(
value=func(bc) if bc.should_extract else bc.value,
bc=bc,
)
mapping.setdefault(fvs, extracted_bc)
return mapping
@cached_property
def reverse_map(self):
r = defaultdict(set)
for k, v in self.map.items():
r[v].add(k)
return r
@cached_property
def fv_set(self):
return set(x[0] for x in self.map.keys())
@cached_property
def offset_partitioning_bc_cells(self):
cells = self.fv_adjacent_cells
cellavg = self.fv_adjacent_cellavg.magnitude.vector()[cells]
return (cells, cellavg)
@cached_property
def fv_adjacent_cells(self):
'''Return list of cells that are adjacent to the facet values.'''
return self.algo_get_fv_adjacent_cells(
facet_function=self.mesh_data.facet_function, fvs=self.fv_set)
@cached_property
def fv_adjacent_cells_function(self):
'''Returns binary cellwise constant (DG0) function, where cell
has value 1 iff cell is in :py:attr:`fv_adjacent_cells`.'''
func = dolfin.Function(self.mesh_util.space.DG0)
func.vector()[:] = 0.0 # TODO: necessary?
func.vector()[self.fv_adjacent_cells] = 1.0
return func
@cached_property
def fv_adjacent_cellavg(self):
'''DG0 function representing the average of the boundary
condition value in adjacent facets.'''
return self.get_fv_adjacent_cellavg()
[docs] def get_fv_adjacent_cellavg(self):
return self.algo_fvs_to_cell_averages(
mesh_util=self.mesh_util,
bcdata=self.reverse_map,
fv_adjacent_cells_function=self.fv_adjacent_cells_function)
[docs] def algo_get_fv_adjacent_cells(self, facet_function, fvs):
'''Get list of cells adjacent to facets that are marked with a
facet value in ``fvs`` inside ``facet_function``.'''
indices = np.nonzero(np.isin(facet_function.array(), list(fvs)))[0]
mesh = facet_function.mesh()
D = mesh.topology().dim()
Facet = dolfin.Facet
cells = []
for facet_index in indices:
facet = Facet(mesh, facet_index)
for cell in facet.entities(D):
cells.append(cell)
return list(sorted(frozenset(cells)))
[docs] def algo_fvs_to_cell_averages(
self, mesh_util, bcdata,
fv_adjacent_cells_function, target_unit=None):
'''Algorithm for turning BC data into cell averages.
Parameters
----------
mesh_util: :py:class:`.mesh_util.MeshUtil`
Mesh utilities and data.
bcdata: dict
Mapping `{expr: {(facet_value, sign), ...}}`.
target_unit: optional
Unit to ``expr`` to in the ``bcdict`` argument.
fv_adjacent_cells_function: :py:class:`dolfin.Function`
DG0 function containing ones and zeros; ones iff the cell is
adjacent to a BC.
'''
mu = mesh_util
dless = mu.unit_registry('dimensionless')
DG0 = mu.space.DG0
dx = mu.dx
u = dolfin.TrialFunction(DG0)
v = dolfin.TestFunction(DG0)
form = DelayedForm()
for expr, fvss in bcdata.items():
if not expr.bc.external:
logger.warning("skipping internal BC {!r}".format(expr.bc))
assert not (expr.bc.internal and expr.bc.should_extract), (
"we don't support default BC value-setting for internal BCs yet"
)
continue
expr = dless * expr.value
if target_unit is None:
target_unit = expr.units
expr_ = expr.m_as(target_unit)
for fv, sign in fvss:
# TODO: figure out what to do with the sign - relevant
# for internal boundary conditions
form += mu.ds(subdomain_id=fv)*v('+')*(u('+') - expr_)
if target_unit is None:
target_unit = dless
form += dx*u*v*(1.0 - fv_adjacent_cells_function)
form = form.delete_units().to_ufl().m
u = dolfin.Function(DG0)
dolfin.solve(dolfin.lhs(form) == dolfin.rhs(form), u, [])
return target_unit*u
[docs] def get_natural_bc(self, ignore_sign=None, bc_reverse_map=None, bc_type='external'):
mu = self.mesh_util
if ignore_sign is None:
ignore_sign = self.natural_ignore_sign
if bc_reverse_map is None:
bc_reverse_map = self.reverse_map
nbc_expr = DelayedForm()
for extracted_bc, fvss in bc_reverse_map.items():
nbc_ds = None
if bc_type=='internal' and extracted_bc.bc.internal:
nbc_ds = mu.dS
elif bc_type=='external' and extracted_bc.bc.external:
nbc_ds = mu.ds
else:
continue
plus_fvs = set()
minus_fvs = set()
for fv, sign in fvss:
if sign > 0:
plus_fvs.add(fv)
else:
minus_fvs.add(fv)
if ignore_sign:
expr = nbc_ds(tuple(sorted(plus_fvs | minus_fvs)))
else:
expr = (nbc_ds(tuple(sorted(plus_fvs))) -
nbc_ds(tuple(sorted(minus_fvs))))
nbc_expr += expr * extracted_bc.value
return nbc_expr
[docs] def get_essential_bc(self, var=None, bc_reverse_map=None):
if var is None:
var = self.essential_var
if bc_reverse_map is None:
bc_reverse_map = self.reverse_map
mu = self.mesh_util
fsr = mu.function_subspace_registry
ebc_space = fsr.get_function_space(
var.m, collapsed=False)
ebc_collapsed_space = fsr.get_function_space(
var.m, collapsed=True)
facet_function = mu.mesh_data.facet_function
var_units = var.units
essential_bcs = []
for value, fvss in bc_reverse_map.items():
value = value.value
value_ = dolfin.project(value.m_as(var_units),
ebc_collapsed_space)
# TODO: determine when usage of geometric marking is necessary
for fv, sign in fvss:
essential_bcs.append(
dolfin.DirichletBC(ebc_space, value_, facet_function, fv))
return essential_bcs
[docs]class Spatial(SetattrInitMixin):
'''Contains spatially dependent quantities, such as material
parameters, as well as boundary conditions.
Parameters
----------
parent: object
Attributes such as mesh_data are taken from here by default.
mesh_util: :py:class:`.mesh_util.MeshUtil`, optional
Mesh-specific utilities. By default taken from parent.
material_to_region_map:
Material to region map. By default taken from
:py:attr:`~mesh_data`.
Attributes
----------
value_rules: defaultdict(SortedList)
Rules for spatially dependent values. Must be dict of list of dicts.
bcs: defaultdict(SortedList)
Boundary conditions.
'''
@property
def mesh_data(self):
return self.mesh_util.mesh_data
@property
def material_to_region_map(self):
return self.mesh_data.material_to_region_map
@property
def mesh_util(self):
return self.parent.mesh_util
@cached_property
def value_rules(self):
return defaultdict(lambda: SortedList(key=lambda x: x.priority))
@cached_property
def bcs(self):
return defaultdict(lambda: SortedList(key=lambda x: x.priority))
[docs] def get(self, key):
'''Get spatially dependent quantity named `key`.'''
mapping = {}
for rule in self.value_rules[key]:
region = rule.region
cvs = self.mesh_data.evaluate_topology(region)
for cv in cvs:
mapping.setdefault(cv, rule.value)
if not len(mapping):
raise AssertionError(
"no value rules for key {!r}, cannot deduce shape or units"
.format(key))
sample_value = rule.value
if hasattr(sample_value, 'units'):
units = sample_value.units
mapping = {k: v.m_as(units) for k, v in mapping.items()}
expr = self.mesh_util.cell_function_function.make_cases(
{k: v for k, v in mapping.items()})
if hasattr(sample_value, 'units'):
expr = units * expr
return expr
[docs] def add_BC(self, key, region, value, priority=0,
require_natural=False, internal=False, external=True, **kwargs):
if isinstance(region, str):
raise AssertionError(
"`region` argument should not be a string, have you "
"swapped the first two arguments by mistake?")
value = constantify_if_literal(value)
bc = BoundaryCondition(region=region,
value=value,
priority=priority,
require_natural=require_natural,
internal=internal,
external=external,
**kwargs)
self.bcs[key].add(bc)
return bc
[docs] def add_rule(self, key, region, value, priority=0):
'''Add a rule for a spatially-dependent value.
Parameters
----------
key: str
Name of the spatial property, e.g., :code:`"temperature"` or
:code:`"poisson/static_rho"`.
region: :py:class:`.topology.CellRegion`
Cell region where this value applies.
value: :py:class:`ufl.core.expr.Expr`
Value.
priority:
Priority of this rule. Lower values have higher priority.
'''
value = constantify_if_literal(value)
self.value_rules[key].add(ValueRule(
region=region,
value=value,
priority=priority))
[docs] def add_value_rule(self, region, key, value, priority=0):
'''Deprecated. Use :py:meth:`add_rule`.'''
warnings.warn(FutureWarning(
"`Spatial.add_value_rule` will be removed soon. Replace it "
"with `Spatial.add_rule`. Note the different argument order!"))
return self.add_rule(
region=region, key=key, value=value, priority=priority)
[docs] def add_material_data(self, key, material_name, value, priority=0):
region = self.material_to_region_map.get(material_name, None)
if region is None:
return
self.add_rule(
key=key,
region=region,
value=value,
priority=priority)
[docs] def make_combined_bc(self, key, value_extractor=None, **kwargs):
cbc = CombinedBoundaryCondition(
bcs=self.bcs[key], value_extractor=value_extractor,
mesh_util=self.mesh_util, **kwargs)
return cbc
[docs] def get_mixed_bc_pair(
self,
essential_key, essential_var,
natural_key , natural_var,
natural_ignore_sign=False,
essential_extract=None,
natural_extract=None):
'''Get mixed BC pair, ensuring that the essential and natural
BCs don't overlap.
Parameters
----------
essential_key: str
Essential BC subfunction key.
natural_key: str
Natural BC subfunction key.
natural_ignore_sign: bool, optional
Whether to ignore facet sign/orientation in the natural
BC. Default is false.
essential_extract: callable, optional
Function that will be called to extract a BC value from a boundary
condition object. Use this to implement any sort of transformation
you want to apply to the boundary condition value.
natural_extract: callable, optional
See `essential_extract` above.
Returns
-------
essential_bcs: list
List of :py:class:`dolfin.DirichletBC`.
natural_bc_expr: UFL expr
Natural BC expression.
'''
if essential_key is not None:
ebc = self.make_combined_bc(
essential_key, essential_extract,
essential_var=essential_var)
else:
ebc = None
if natural_key is not None:
nbc = self.make_combined_bc(
natural_key, natural_extract,
natural_ignore_sign=natural_ignore_sign)
else:
nbc = None
if ebc and nbc and (ebc.fv_set & nbc.fv_set):
raise AssertionError(
"Essential and natural boundary conditions overlap. "
"Keys {!r} and {!r}.".format(essential_key, natural_key))
return (ebc, nbc)
[docs] def get_single_essential_bc(
self, key, var, extract=None):
'''Get single essential BC.
Parameters
----------
key: str
Essential BC key.
var: str
Essential BC variable.
extract: callable, optional
See :code:`essential_extract` under :py:meth:`get_mixed_bc_pair`.
Returns
-------
essential_bcs: list
List of :py:class:`dolfin.DirichletBC`.
'''
return self.get_mixed_bc_pair(
essential_key=key,
essential_var=var,
natural_key=None,
natural_var=None,
essential_extract=extract)[0]
[docs] def get_vanilla_bc_pair(self, variable_key):
'''not implemented because all our methods are mixed methods'''
raise NotImplementedError()