It is important to realize that a linear PDE with nonlinear BCs gives rise to a nonlinear problem, the linearity of the PDE notwithstanding.

Writing Neumann BCs in Galerkin formulations is extremely simple. After integration by parts, we will have boundary terms involving test functions times derivatives of our unknowns. To apply a Neumann BC, we simply replace the derivatives appearing in those boundary terms with their value determined by the BC. Homogeneous Neumann BCs are a particularly simple special case because the boundary integrals go to zero, letting us satisfy the BC by simply *ignoring* the boundary integrals on those regions! It doesn't get much simpler than that.

In many cases we can enforce Dirichlet boundary conditions on certain nodes or edges by replacing the PDE at those entities with an equation representing the Dirichlet BCs. This can be done nodewise or weakly over the boundary cell; note that nodewise BC application is a special case of weak BC application obtained by using a nodal quadrature rule. A justification of this scheme is given below in section Optional: More about Dirichlet BCs.

In Sundance we write such a boundary condition using the EssentialBC object. The most general constructor for an EssentialBC takes as arguments

- the CellFilter on which the BC is to be applied
- a weak form of the BC.
- a QuadratureFamily

On the boundary regions , , and we have boundary conditions

We form a weak equation by multiplying by a test function and integrating, giving

which must hold for all in a suitable subspace. On the Neumann and Robin regions, we can substitute and , respectively, for the normal derivative , so that these BCs are simply included in the weak form as written. The method for applying Dirichlet BCs is described in Optional: More about Dirichlet BCs.

This equation is written in Sundance with code such as the following

// Define test and unknown functions Expr v = new TestFunction(some basis); Expr u = new UnknownFunction(some basis); // Define a cell filter for the interior of the region CellFilter omega = new MaximalCellFilter(); // Define cell filters for the surfaces on which // neumann, robin, and dirichlet are to be applied. CellFilter gamma_neumann = [some cell filter definition]; CellFilter gamma_robin = [some other cell filter definition]; CellFilter gamma_diri = [yet another cell filter definition]; // Define an expression for the flux on the Neumann surface Expr G = [some expression]; // Define expressions for the variables in the Robin BC Expr alpha = [some expression]; Expr u_R = [some other expression]; // Define expressions for the variable in the Dirichlet BC Expr u_D = [some expression]; Expr eqn = Integral(omega, (grad*v)*(grad*u)) + Integral(gamma_neumann, v*G) + Integral(gamma_robin, v*alpha*(u-u_R)); Expr essBC = EssentialBC(gamma_diri, v*(u-u_D));

Without loss of generality, we'll use as a model problem the Laplace equation with Dirichlet conditions on the entire boundary:

As usual, we form a weak equation by multiplying by a test function and integrating, giving

which must hold for all in a suitable subspace. We now impose a Robin boundary condition

which in the limit approaches the Dirichlet condition . With this BC, the weak equation becomes

At this point, we could simply take a small but nonzero and solve the problem. This is equivalent to a penalty method for imposing the BC, and unfortunately leads to poorly conditioned linear systems. We want to do better, and take the limit so that the Dirichlet BC is satisfied exactly, yet we also want to avoid ill-conditioning. As we will show, we can often have it both ways, taking the limit without ill-conditioning, but the argument must be developed carefully because the way in which that limit is taken is important.

We will discretize the system before taking the limit. We expand in a series of basis functions,

and evaluate Blah} with each member of as test functions. The discrete equations are

where the matrices are defined as

At this point we make a critical assumption: that we are using either Lagrange or Serendipity basis functions. These basis functions are { nodal}, by which we mean we can make an association between basis functions and mesh nodes, and then define the basis functions in such a way that the -th basis function is zero at every node but the -th. Furthermore, these basis sets have the important property that { only those basis functions associated with nodes on the surface have nonzero values on the surface .} We can now partition the equations and variables according to whether each row or column is associated with a node on . Notice that the dangerous denominator appears only in terms involving , whose elements are obtained by integrating over . Thus all elements of involving nodes not on are zero. Marking elements associated with with the superscript and all others with the superscript , we can write each matrix as a block matrix, e.g.,

In particular, from the argument above we know that all off-boundary elements of are zero,

With this notation, we can write out the system,

Multiplying the second row by $$ removes the explosive denominator

so that it is safe to take the limit $ 0$, yielding the well-behaved system

This entire process can be summarized as such: for surfaces on which Dirichlet BCs apply, we simply * replace* all equations by the Dirichlet boundary term. The bookkeeping for this would be daunting if we had to keep track of off-boundary and on-boundary nodes in the equation specification. However, in Sundance the specification that a term is to be handled in this special way via replacement (or more precisely via limit and rescaling) is done by putting the term inside an EssentialBC object rather than an Integral object. The Sundance computational kernel will figure out which rows get replaced and rescaled.

- Pick one of the two boundary conditions and apply it at the corner point.
- Average the two boundary conditions.

It is often possible to find a local approximation to a nonlocal boundary condition, which will allow you to proceed at the prce of introducing some modeling error.