Solving Poisson equation with current BC using FEM

Question

I have a 3-dimensional domain D, with 3 types of BC which I am trying to solve the Poisson equation on. $$-\nabla \cdot (\sigma(x,y,z)\nabla u)=0$$
Insulator (Neumann BC)
Electrode set at some potential (Dirichlet BC).
Additionally I have some electrodes $E_l$ which are set to have a fixed current $I_l$, so on these the boundary condition is $\int_{E_l} \sigma(x,y,z)(du/dn) dA=I_l$.

This boundary condition can given a mesh be translated into equations of the form $\sum a_iu_i=0$, where $u_i$ is the potential at node i and $a_i$ are known mesh dependend constants. But I dont know how to add this equation to the matrix to be solved, or what equations to replace it with.

I cannot assume that the current density is constant on each electrode, for instance some electrodes have net current 0 but current flowing in on one side and out the other, I know commercial software solves this problem just as efficiently as with regular BC.

How would I go about incoprorating this third type of BC into FEM? Any description or references is what im looking for.

Answer 1

It is standard procedure for general case like that. You can enforce condition for electrode by Lagrange multipliers, $$ L(u,\lambda) = \int_\Omega \sigma u_{,i} u_{,i} \, \textrm{d}V - \int_{\partial\Omega_\sigma} u n_i j_i \textrm{d}S + \lambda\left(\int_E \sigma u_{,n}\, \textrm{d}S - I \right) $$ that for discretised problem is $$ \left[ \begin{array}{cc} A & C^T\\ C & 0 \end{array} \right] \left\{ \begin{array}{c} u \\ \lambda \end{array} \right\} =\left[ \begin{array}{c} f \\ I \end{array} \right] $$ where $$ C_i = \int_E \sigma n_j \phi_{i,j} \textrm{d}S $$ and $\phi_i$ is shape functions at the node. The technical difficulty here is to solve for the gradient of the shape functions on the surface, but that's solvable with good finite element library. You have to enforce Dirichlet bc a priori in standard way.

Note if you have more electrodes, that comes with more Lagrange multipliers. You can efficiently solve the problem using Schur complement and block solvers like http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/PC/PCFIELDSPLIT.html

EDIT:

Another way is to use projection matrices, avoiding explicitly Lagrange multipliers, once you have constrains matrix $Cu=g$, you can solve problem $$ C^TC + Q^TKQ = C^Tg + Q^T(f-KRg) $$ where $$ Q = I-P^TKP $$, $$ P = C^T(CC^T)^{-1}C $$ $$ R = (CC^T)^{-1}C $$

You can find details of this in

Mark Ainsworth. Essential boundary conditions and multi-point constraints in finite element analysis. Computer Methods in Applied Mechanics and Engineering, 190(48):6323–6339, 2001.

I have been using this couple of times, you can implement this such that you do not have explicitly $Q^TKQ$, you can see implementation here.

Answer 2

Like @likask points out in the other answer, using Lagrange multipliers is one way to do this. The other is to recognize that if you have constraints of the form $\sum_i a_i U_i=0$, then you can select one degree of freedom (typically the one with the largest $a_i$, to maintain numerical stability and rewrite the constraint as $$ U_\ast = -\sum_{i, i\neq \ast} \frac{a_i}{a_\ast} U_i, $$ i.e., one degree of freedom is a linear combination of a bunch of others. This is then no different from what you do, for example, with hanging nodes, with tangential flow boundary conditions, and a bunch of other things in finite element codes. There is a well-developed machinery to deal with such constraints.

An example where this is done in practice is the step-11 tutorial program of deal.II. If you're unsure how constraints such as this are handled algorithmically, you may be interested in my video lecture 16.