Skip to content

Commit

Permalink
Merge pull request #2007 from firedrakeproject/wence/feature/no-geom-bcs
Browse files Browse the repository at this point in the history 
Interior bcs
  • Loading branch information
@wence-
wence- authored May 13, 2021
2 parents eada709 + fc6bc46 commit 31b3c11
Show file tree
Hide file tree
Showing 23 changed files with 1,165 additions and 362 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,3 +38,4 @@ __pycache__/
/firedrake/cython/supermeshimpl.c
/firedrake_configuration/configuration.json
/firedrake.egg-info
/docs/source/element_list.csv
59 changes: 33 additions & 26 deletions docs/source/variational-problems.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -520,6 +520,34 @@ details of how Firedrake applies strong boundary conditions are
slightly involved and therefore have :doc:`their own section
<boundary_conditions>` in the manual.

Boundary conditions on interior facets
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

If you wish to apply strong boundary conditions to interior facets of
your mesh, this is transparently supported. You should arrange that
your mesh generator marks those facets on which you wish to apply
boundary conditions, and just use the subdomain ids as usual.

Special subdomain ids
~~~~~~~~~~~~~~~~~~~~~

As well as integer subdomain ids that come from marked portions of the
mesh, Firedrake also supports the magic string ``"on_boundary"`` to
apply a boundary condition to all exterior facets of the mesh.
Further, on :doc`:extruded meshes <extruded-meshes>` the special
strings ``"top"`` and ``"bottom"`` can be used to apply a boundary
condition on respectively the top and bottom of the extruded domain.

.. note::

These special strings cannot be combined with integer ids, so if
you want to apply boundary data on an extruded mesh on (say) ids
``1`` and ``2`` as well as the top of the domain you would write

.. code-block:: python3
bcs = [DirichletBC(V, ..., (1, 2)), DirichletBC(V, ..., "top")]
Specifying conditions on components of a space
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Expand Down Expand Up @@ -559,35 +587,14 @@ velocity is zero, using the same function space definitions:
bcv_x = DirichletBC(W.sub(0).sub(0), Constant(0), boundary_ids)
.. note::

Extruded meshes have full support for indexing
:py:class:`~.MixedFunctionSpace`\s, but currently do not support
indexing on :py:func:`~.VectorFunctionSpace`\s.

Boundary conditions in discontinuous spaces
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

The default method Firedrake uses to determine where to apply strong
boundary conditions is :py:data:`"topological"`, meaning that nodes
topologically associated with a boundary facet will be included. In
discontinuous spaces, however, the nodes to be included do not all
live on boundary facets, in this case, you should use the
:py:data:`"geometric"` method for determining boundary condition
nodes. In this case, nodes associated with basis functions that do
not vanish on the boundary are included. This method can be used to
impose strong boundary conditions on discontinuous galerkin spaces, or
no-slip conditions on :math:`H(\textrm{div})` spaces. To select the method used for
determining boundary condition nodes, use the :py:data:`method`
argument to the :py:class:`DirichletBC` constructor. For example, to
select geometric boundary node determination we would write:

.. code-block:: python3
V = FunctionSpace(mesh, 'DG', 2)
bc = DirichletBC(V, 1.0, subdomain_id, method="geometric")
...
Firedrake uses the topological association of nodes to facets to
determine where to apply strong boundary conditions. For spaces where
nodes are not topologically associated with the boundary facets, such
as discontinuous Galerkin spaces, you should instead apply boundary
conditions weakly.

Time dependent boundary conditions
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
Expand Down
72 changes: 37 additions & 35 deletions firedrake/bcs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,14 +35,8 @@ class BCBase(object):
to indicate all of the boundaries of the domain. In the case of extrusion
the ``top`` and ``bottom`` strings are used to flag the bcs application on
the top and bottom boundaries of the extruded mesh respectively.
:arg method: the method for determining boundary nodes. The default is
"topological", indicating that nodes topologically associated with a
boundary facet will be included. The alternative value is "geometric",
which indicates that nodes associated with basis functions which do not
vanish on the boundary will be included. This can be used to impose
strong boundary conditions on DG spaces, or no-slip conditions on HDiv spaces.
'''
def __init__(self, V, sub_domain, method="topological"):
def __init__(self, V, sub_domain):

# First, we bail out on zany elements. We don't know how to do BC's for them.
if isinstance(V.finat_element, (finat.Argyris, finat.Morley, finat.Bell)) or \
Expand All @@ -51,9 +45,6 @@ def __init__(self, V, sub_domain, method="topological"):
self._function_space = V
self.comm = V.comm
self.sub_domain = sub_domain
if method not in ["topological", "geometric"]:
raise ValueError("Unknown boundary condition method %s" % method)
self.method = method
# If this BC is defined on a subspace (IndexedFunctionSpace or
# ComponentFunctionSpace, possibly recursively), pull out the appropriate
# indices.
Expand Down Expand Up @@ -151,7 +142,7 @@ def hermite_stride(bcnodes):
bcnodes = []
for s in sub_d:
if isinstance(s, str):
bcnodes.append(hermite_stride(self._function_space.boundary_nodes(s, self.method)))
bcnodes.append(hermite_stride(self._function_space.boundary_nodes(s)))
else:
# s is of one of the following formats:
# facet: (i, )
Expand All @@ -166,7 +157,7 @@ def hermite_stride(bcnodes):
# intersection of facets
# Edge conditions have only been tested with Lagrange elements.
# Need to expand the list.
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss, self.method)))
bcnodes1.append(hermite_stride(self._function_space.boundary_nodes(ss)))
bcnodes1 = functools.reduce(np.intersect1d, bcnodes1)
bcnodes.append(bcnodes1)
return np.concatenate(bcnodes)
Expand Down Expand Up @@ -243,6 +234,13 @@ def extract_forms(self, form_type):
class DirichletBC(BCBase, DirichletBCMixin):
r'''Implementation of a strong Dirichlet boundary condition.
.. note:
This uses facet markers in the domain, so may be used to
applied strong boundary conditions to interior facets (if they
have an appropriate mesh marker). The "on_boundary" string only
applies to the exterior boundaries of the domain.
:arg V: the :class:`.FunctionSpace` on which the boundary condition
should be applied.
:arg g: the boundary condition values. This can be a :class:`.Function` on
Expand All @@ -257,17 +255,25 @@ class DirichletBC(BCBase, DirichletBCMixin):
to indicate all of the boundaries of the domain. In the case of extrusion
the ``top`` and ``bottom`` strings are used to flag the bcs application on
the top and bottom boundaries of the extruded mesh respectively.
:arg method: the method for determining boundary nodes. The default is
"topological", indicating that nodes topologically associated with a
boundary facet will be included. The alternative value is "geometric",
which indicates that nodes associated with basis functions which do not
vanish on the boundary will be included. This can be used to impose
strong boundary conditions on DG spaces, or no-slip conditions on HDiv spaces.
:arg method: the method for determining boundary nodes.
DEPRECATED. The only way boundary nodes are identified is by
topological association.
'''

@DirichletBCMixin._ad_annotate_init
def __init__(self, V, g, sub_domain, method="topological"):
super().__init__(V, sub_domain, method=method)
def __init__(self, V, g, sub_domain, method=None):
if method == "geometric":
raise NotImplementedError("'geometric' bcs are no longer implemented. Please enforce them weakly")
if method not in {None, "topological"}:
raise ValueError(f"Unhandled boundary condition method '{method}'")
if method is not None:
import warnings
with warnings.catch_warnings():
warnings.simplefilter('always', DeprecationWarning)
warnings.warn("Selecting a bcs method is deprecated. Only topological association is supported",
DeprecationWarning)
super().__init__(V, sub_domain)
if len(V) > 1:
raise ValueError("Cannot apply boundary conditions on mixed spaces directly.\n"
"Apply to the components by indexing the space with .sub(...)")
Expand All @@ -286,16 +292,14 @@ def function_arg(self):
self._function_arg_update()
return self._function_arg

def reconstruct(self, field=None, V=None, g=None, sub_domain=None, method=None, use_split=False):
def reconstruct(self, field=None, V=None, g=None, sub_domain=None, use_split=False):
fs = self.function_space()
if V is None:
V = fs
if g is None:
g = self._original_arg
if sub_domain is None:
sub_domain = self.sub_domain
if method is None:
method = self.method
if field is not None:
assert V is not None, "`V` can not be `None` when `field` is not `None`"
V = self.as_subspace(field, V, use_split)
Expand All @@ -307,9 +311,9 @@ def reconstruct(self, field=None, V=None, g=None, sub_domain=None, method=None,
(V.parent is None or V.parent.parent == fs.parent.parent) and \
(V.parent is None or V.parent.index == fs.parent.index) and \
g == self._original_arg and \
sub_domain == self.sub_domain and method == self.method:
sub_domain == self.sub_domain:
return self
return type(self)(V, g, sub_domain, method=method)
return type(self)(V, g, sub_domain)

@function_arg.setter
def function_arg(self, g):
Expand Down Expand Up @@ -441,14 +445,13 @@ class EquationBC(object):
:param Jp: a form used for preconditioning the linear system,
optional, if not supplied then the Jacobian itself
will be used.
:arg method: see :class:`.DirichletBC` (optional)
:arg V: the :class:`.FunctionSpace` on which
the equation boundary condition is applied (optional)
:arg is_linear: this flag is used only with the `reconstruct` method
:arg Jp_eq_J: this flag is used only with the `reconstruct` method
'''

def __init__(self, *args, bcs=None, J=None, Jp=None, method="topological", V=None, is_linear=False, Jp_eq_J=False):
def __init__(self, *args, bcs=None, J=None, Jp=None, V=None, is_linear=False, Jp_eq_J=False):
from firedrake.variational_solver import check_pde_args, is_form_consistent
if isinstance(args[0], ufl.classes.Equation):
# initial construction from equation
Expand Down Expand Up @@ -487,9 +490,9 @@ def __init__(self, *args, bcs=None, J=None, Jp=None, method="topological", V=Non
# Argument checking
check_pde_args(F, J, Jp)
# EquationBCSplit objects for `F`, `J`, and `Jp`
self._F = EquationBCSplit(F, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._F for bc in bcs], method=method, V=V)
self._J = EquationBCSplit(J, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._J for bc in bcs], method=method, V=V)
self._Jp = EquationBCSplit(Jp, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._Jp for bc in bcs], method=method, V=V)
self._F = EquationBCSplit(F, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._F for bc in bcs], V=V)
self._J = EquationBCSplit(J, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._J for bc in bcs], V=V)
self._Jp = EquationBCSplit(Jp, u, sub_domain, bcs=[bc if isinstance(bc, DirichletBC) else bc._Jp for bc in bcs], V=V)
elif all(isinstance(args[i], EquationBCSplit) for i in range(3)):
# reconstruction for splitting `solving_utils.split`
self.Jp_eq_J = Jp_eq_J
Expand Down Expand Up @@ -533,12 +536,11 @@ class EquationBCSplit(BCBase):
:arg sub_domain: see :class:`.DirichletBC`.
:arg bcs: a list of :class:`.DirichletBC`s and/or :class:`.EquationBC`s
to be applied to this boundary condition equation (optional)
:arg method: see :class:`.DirichletBC` (optional)
:arg V: the :class:`.FunctionSpace` on which
the equation boundary condition is applied (optional)
'''

def __init__(self, form, u, sub_domain, bcs=None, method="topological", V=None):
def __init__(self, form, u, sub_domain, bcs=None, V=None):
# This nested structure will enable recursive application of boundary conditions.
#
# def _assemble(..., bcs, ...)
Expand All @@ -551,7 +553,7 @@ def __init__(self, form, u, sub_domain, bcs=None, method="topological", V=None):
self.u = u
if V is None:
V = self.f.arguments()[0].function_space()
super(EquationBCSplit, self).__init__(V, sub_domain, method=method)
super(EquationBCSplit, self).__init__(V, sub_domain)
# overwrite bcs
self.bcs = bcs or []
for bc in self.bcs:
Expand Down Expand Up @@ -602,10 +604,10 @@ def reconstruct(self, field=None, V=None, subu=None, u=None, row_field=None, col
if action_x is not None:
assert len(form.arguments()) == 2, "rank of self.f must be 2 when using action_x parameter"
form = ufl_expr.action(form, action_x)
ebc = EquationBCSplit(form, subu, self.sub_domain, method=self.method, V=W)
ebc = EquationBCSplit(form, subu, self.sub_domain, V=W)
for bc in self.bcs:
if isinstance(bc, DirichletBC):
ebc.add(bc.reconstruct(V=W, g=bc.function_arg, sub_domain=bc.sub_domain, method=bc.method, use_split=use_split))
ebc.add(bc.reconstruct(V=W, g=bc.function_arg, sub_domain=bc.sub_domain, use_split=use_split))
elif isinstance(bc, EquationBCSplit):
bc_temp = bc.reconstruct(field=field, V=V, subu=subu, u=u, row_field=row_field, col_field=col_field, action_x=action_x, use_split=use_split)
# Due to the "if index", bc_temp can be None
Expand Down
Loading

0 comments on commit 31b3c11

Please sign in to comment.