.. |fenics| replace:: FEniCS .. _basic-bcs: How to define and apply boundary conditions? ============================================ .. note:: This `tutorial `_ covers the application of different kind of boundary conditions (Dirichlet, Neumann and Robin) following different strategies (from the basic use of functions to define boundaries, to more complex approaches as using compiled subdomains). Please refer to it for a more detailed overview. By definition, Dirichlet boundary conditions represent degrees of freedom (dofs) for which we already know the solution. Therefore, we essentially need to provide |fenics| with the corresponding dofs or a way to find the corresponding dofs (e.g. a function that defines if a point belongs to the Dirichlet boundary), and the corresponding values. Homogeneous Neumann boundary conditions are naturally handled by the finite element method (meaning boundaries without explicitly defined boundary conditions are considered homogeneous Neumann). Inhomogeneous Neumann boundary conditions are represented by an additional term in the variational approach (which will be converted in an additional vector belonging to the righ-hand side of the variational problem, i.e. to the linear form). .. note:: In the strategy explained above, Dirichlet boundary conditions are strongly imposed. An alternative is to impose them `weakly `_. Dirichlet boundary conditions ----------------------------- .. _define-bc: Define (a simple example) ************************* Imagine we want to apply a Dirichlet boundary condition to the left edge of a unit square. First, we need to define a function that identifies if a given point belongs to the edge: .. code-block:: python from dolfin.cpp.math import near def boundary(x, on_boundary): return on_boundary and near(x[0], 0.) ``x`` is an array with the coordinates of the point. ``on_boundary`` is a boolean that is ``True`` if the point belongs to the boundary. ``near`` compares two floats (by default a ``tol= 3e-16``). Then, we create a Dirichlet boundary condition object: .. code-block:: python from dolfin.fem.dirichletbc import DirichletBC bc = DirichletBC(V, u_D, boundary) ``V`` is the function space and ``u_D`` the desired value. More complex boundaries are hard to define using simple mathematical expressions. Boundary markers are then the simplest way to identify them (check :ref:`this ` out for more information). Apply: black-box case ********************** The black-box approach to apply Dirichlet boundary conditions consists in simply passing them to the ``solve`` function: .. code-block:: python solve(a == L, u, bc, solver_parameters=solver_parameters) If you are solving a symmetric problem, ensure to pass ``symmetric = True`` in ``solver_parameters``, in order to keep symmetry after the application of boundary conditions and use a convenient solver (more information on how it impacts solver choice :ref:`here `). Apply: breaking symmetry ************************ On the other hand, it is usually useful to apply boundary conditions directly on assembled matrices. The simplest (yet not recommended) way to do it consists in zeroing-out rows that represent Dirichlet dofs. That's the task performed by ``bc.apply``. As example, let's take the bilinear form of Poisson equation: .. math:: a(u, v) = \int_{\Omega} \nabla u \cdot \nabla v \mathrm{~d} x, which is translated in ``ufl`` language as (``V`` is the function space): .. code-block:: python from dolfin.function.argument import TrialFunction from dolfin.function.argument import TestFunction from ufl import dot from ufl import grad from ufl import dx u_h = TrialFunction(V) v_h = TestFunction(V) a = dot(grad(u_h), grad(v_h))*dx To assemble the matrix: .. code-block:: python from dolfin.fem.assembling import assemble A = assemble(a) Finally, to apply boundary conditions: .. code-block:: python for bc in bcs: bc.apply(A) A pictorial representation of the assembled matrices before (left) and after (right) boundary condition application is shown next (yellow represents non-zero elements). .. image:: images/assembled_matrix.png :width: 45% .. image:: images/assembled_matrix_bcs.png :width: 45% | The drawback of this approach is that initially symmetric matrices become non-symmetric (e.g. take a look to first row and column)... and symmetry is a property we want to keep at all costs! In the right-hand side, this solution consists in replacing the components corresponding to Dirichlet dofs by the desired values (the Dirichelet dofs equations become :math:`u_i = u_{D_i}`). Apply: keeping symmetry *********************** A better way to apply boundary conditions is to use ``assemble_system``, since it keeps symmetry when applying boundary conditions. For that, we need both the bilinear and linear forms of the problem. For Poisson equation, the linear form is: .. math:: L(v) = \int_{\Omega} f v \mathrm{~d} x, which is translated in ``ufl`` language as (``V`` is the function space): .. code-block:: python from dolfin.function.constant import Constant from ufl import dx f = Constant(-6.) # just as example L = f*v_h*dx To assemble the matrix and the vector: .. code-block:: python from dolfin.fem.assembling import assemble_system A, b = assemble_system(a, L, bcs) This way, ``A`` is kept symmetric, as the next pictorial representation shows. .. image:: images/assembled_matrix_sym.png :align: center :width: 45% | ``assemble_system`` incorporates (in a symmetric way) the boundary conditions directly in the element matrices and vectors, prior to assembly (check out `Logg's book `_ to learn more about the assembly algorithm). The advantage of taking care of boundary conditions at element level is that the process of eliminating whole rows and columns is more `efficient `_ (than when matrices are stored in compressed row format). Apply: keeping symmetry but varying right-hand side (advanced) ************************************************************** For time-dependent problems, it is usually convenient to assemble the left-hand side only once (as the assembly operation can be costly for large problems). Following the "breaking symmetry" approach, it is trivial to achieve this goal, as the application of boundary conditions is done independently for each equation side. For the "keeping symmetric" approach, things get slightly more complicated: it is necessary to keep track of an additional vector that stores the results of the elimination process. The process of assembling the left-hand side is: .. code-block:: python dummy = Constant(0.)*v_h*dx # integrates to zero A, dummy_vec = assemble_system(a, dummy) Mathematically, this dummy vector is :math:`\sum_i u_{D_i} a_i`, where :math:`u_{D_i}` represents the value of the Dirichlet boundary condition at dof :math:`i` and :math:`a_i` represents the :math:`i` column of matrix :math:`A` after zeroing-out all the Dirichlet dofs rows. To obtain the right-hand side: .. code-block:: python b_tmp = assemble(L) b = b_tmp + dummy_vec for bc in bcs: bc.apply(b) ``bc.apply(b)`` is performed to ensure the correct values at Dirichlet dofs components. .. note:: Althought quite general, this example ignores time-dependent Dirichlet boundary conditions. Can you see what needs to change? (Tip: writing :math:`u_D` as :math:`u_D(x, t) = u_D(t) u_D(x) + u_D(t) + u_D(x)` can help)