# Implementing a potential problem with a current condition

**URL:** <https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491>\
**Category:** variational formulation\
**Created:** [April 30, 2024, 3:35pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491 "2024-04-30T15:35:41Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [April 30, 2024, 3:35pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/1 "2024-04-30T15:35:41Z")

</div>

![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/f/f64a91d0045c43731e6250e0ec617ca3ceebbf5c.png)

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [April 30, 2024, 4:03pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/2 "2024-04-30T16:03:32Z")

</div>

This doesn’t seem to me to be a well posed problem. The solution u is overconstrained on Gamma\_2, because it is asked to be both equal to V (a Dirichlet condition, which would be acceptable on its own), and such that the average of a partial derivative is equal to one.

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [April 30, 2024, 4:18pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/3 "2024-04-30T16:18:30Z")

</div>

The point here is that V is an unknown quantity, which must be solved for using the Integral constraint involving the gamma term.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [April 30, 2024, 4:48pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/4 "2024-04-30T16:48:34Z")

</div>

OK, I didn’t get that from reading the text the first time. I don’t have specific experience for this kind of constraint, but you may want to start by searching the forum with “nullspace” as keyword, and start reading recent posts about that.

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [April 30, 2024, 6:32pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/5 "2024-04-30T18:32:37Z")

</div>

I found the following post that seems to solve a similar type of problem  
[https://fenicsproject.org/olddocs/dolfin/2018.1.0/python/demos/singular-poisson/demo\_singular-poisson.py.html](https://fenicsproject.org/olddocs/dolfin/2018.1.0/python/demos/singular-poisson/demo_singular-poisson.py.html)  
Unfortunately, the solution is in fenics instead of fenicsx, so I don’t know how to convert it. (I am a complete beginner and trying to teach myself this stuff.) So here is a comment I would like to make about the problem that is outlined, called the “singular Poisson problem”. If one were to choose any single point on the boundary and impose a Dirichlet condition at that point, the problem would no longer be singular. Now impose the condition that the average of u over the entire domain must be zero, and use that condition to solve for the unknown value in the Dirichlet condition. This is completely analogous to what I am trying to do in my problem. Apparently this can be handled in fenics, but how do we do it in fenicsx?

In the tutorial, they use something called a Krylov solver to impose the average condition. My method above seems more straightforward. The question is whether or not fenicsx is set up to handle my approach.

I just noticed that the singular Poisson problem was implemented in dolfinx in

> [@Issues Solving Pure Neumann Problems in dolfinx](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/3):
>
> Thanks @dokken. I was able to replicate the Singular Poisson problem in dolfinx using your help. The dolfin example is [here.](https://fenicsproject.org/olddocs/dolfin/2018.1.0/python/demos/singular-poisson/demo_singular-poisson.py.html) It might be nice to include this example in the dolfinx tutorial. Either way, I’ll just post the code I was able to get working here (apologies if this already exists elsewhere, it just took me a bit to put the pieces together): from mpi4py import MPI from dolfinx import mesh,plot from dolfinx.fem import FunctionSpace,Function,form from dolfinx.fem.petsc import create\_v…

Not sure if I can figure out how to modify this approach to solve the above problem with a current condition, but I will spend some time thinking about it and report my questions or answers.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [May 1, 2024, 12:55pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/6 "2024-05-01T12:55:51Z")

</div>

Having seen your recent posts in [Issues Solving Pure Neumann Problems in dolfinx - #7 by Daniel.Baker](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/7), I’d like to remark that one difference between the problem here and the one in the other post is that here you actually have two constraints:

- the integral constraint on Gamma\_2
- the fact that the solution u must have the same value on all DOFs on Gamma\_2

Having said that, I have no experience in how to include these constraints in dolfinx, and hence I can’t help any further.

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 1, 2024, 3:44pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/7 "2024-05-01T15:44:34Z")

</div>

I agree entirely with your observation. Having said that, it is perhaps worth being more explicit about the second constraint. Let’s choose an arbitrary point x\_0 in Gamma\_2. Then the second constraint says that the difference u(x)-u(x\_0)=0 for every point x in Gamma\_2. Thus, if there are N DOFs for u in Gamma\_2, there will be N-1 constraints of the second type and the integral constraint is also a single constraint. So the total number of constraints is N, the number of DOFs for u in Gamma2. Thus the problem is well defined.

In [Issues Solving Pure Neumann Problems in dolfinx](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/7) one can also choose an arbitrary point x\_0 on the boundary of the region and stipulate a Dirichlet condition at x\_0. Everywhere else is Neumann boundary conditions. The value of the Dirichlet condition is then determined by the integral constraint on u over the entire domain. In this regard, the two problems are very similar.

---

<div class="post-metadata">

**Author:** ![francesco-ballarin](https://avatars.discourse-cdn.com/v4/letter/f/f05b48/32.png) [@francesco-ballarin](https://fenicsproject.discourse.group/u/francesco-ballarin)\
**Post date:** [May 1, 2024, 4:19pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/8 "2024-05-01T16:19:39Z")

</div>

@Daniel.Baker and @dokken may want to continue here the discussion from [Issues Solving Pure Neumann Problems in dolfinx - #8 by dokken](https://fenicsproject.discourse.group/t/issues-solving-pure-neumann-problems-in-dolfinx/10468/8)

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 1, 2024, 8:20pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/9 "2024-05-01T20:20:50Z")

</div>

I will try to reply to comments from @dokken that seem to be relevant to the problem I posed at the top of this thread. I have looked through his suggestions, and I think that, in principle, the information he provided could be used to solve the above problem, if I only understood how to implement it in code. So here is the reference from @dokken that I am interested in

> **[GitHub - jorgensd/dolfinx\_mpc: Extension for dolfinx to handle multi-point...](https://github.com/jorgensd/dolfinx_mpc)**
>
> Extension for dolfinx to handle multi-point constraints. - jorgensd/dolfinx\_mpc

Specifically, he refers to

 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/b/b174c78218ae3e41889a93198719d6d8312d34c8.png)  
In my case, let u\_0, u\_1, …,u\_N be the variables on the boundary Gamma2. These constraints would be  
u\_i=u\_0 for all 0\<i\<N  
There is one more inhomogeneous constraint I need, which, in variational form, could be written as  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/6/6749af8338d8db24358e5f8c89cc5847c69a8c14.png)  
So I think that all of this would fit into the format that it described. My job is to understand how this code is implemented, for example, in the tutorial  
[https://jsdokken.com/dolfinx\_mpc/python/demos/demo\_stokes.html](https://jsdokken.com/dolfinx_mpc/python/demos/demo_stokes.html)  
It would greatly help me to understand what is in this tutorial, if someone could document for me the mathematical PDE’s and boundary conditions that are being solved in it. In the tutorial itself, all I see is code, which so far I have been unable to understand. Thanks ahead of time for your help!

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 2, 2024, 5:12pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/10 "2024-05-02T17:12:09Z")

</div>

This was my bad. Full documentation of the equations is cited in the demo\_stokes in the link  
[https://arxiv.org/pdf/2001.10639.pdf](https://arxiv.org/pdf/2001.10639.pdf)  
It might help to put this reference at the top of the demo. The following reference is also helpful  
[https://hal.science/hal-04077986/document](https://hal.science/hal-04077986/document)

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 2, 2024, 6:24pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/11 "2024-05-02T18:24:30Z")

</div>

A presentation of DOLFINx\_mpc as a framework was done at FEniCS 21: [https://fenics2021.com/slides/dokken.pdf](https://fenics2021.com/slides/dokken.pdf)  
and at Sorbonne/D’Alembert institute fall of 2023: [https://www.youtube.com/watch?v=4LSWOdyiGH8&ab\_channel=InstitutJeanLeRondd’Alembert](https://www.youtube.com/watch?v=4LSWOdyiGH8&ab_channel=InstitutJeanLeRondd%27Alembert)

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 2, 2024, 8:58pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/12 "2024-05-02T20:58:10Z")

</div>

I have dolfinx up and running. How do I find dolfinx\_mpc so I can install it? Sorry to have to keep asking these questions. I should be starting to know this stuff by now.

Got it! I’m up and running.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 2, 2024, 9:06pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/13 "2024-05-02T21:06:20Z")

</div>

Depending on how you installed dolfinx, choose the appropriate way:

> **[GitHub - jorgensd/dolfinx\_mpc: Extension for dolfinx to handle multi-point...](https://github.com/jorgensd/dolfinx_mpc?tab=readme-ov-file#installation)**
>
> Extension for dolfinx to handle multi-point constraints. - jorgensd/dolfinx\_mpc

or

> **[Debian -- Details of package python3-dolfinx-mpc in sid](https://packages.debian.org/sid/python3-dolfinx-mpc)**
>
> Multi-point constraints with FEniCS-X

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 5, 2024, 1:37pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/14 "2024-05-05T13:37:28Z")

</div>

I have some extra questions about dolfinx\_mpc.

1. Will it handle inhomogeneous constraints? For example, an inhomogeneous periodic boundary condition, where the difference between the left and right boundary is a known function?
2. If the known function in question (1) depends on a single unknown constant, can one introduce an additional constraint to determine the value of this constant? For example if the two boundaries are at x=0 and x=1, suppose we want to introduce the condition  
u(x=1)=u(x=0)+K  
where K is to be determined by \int\_Omega u dx =1 (again, an inhomogeneous constraint)?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 5, 2024, 2:55pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/15 "2024-05-05T14:55:13Z")

</div>

No, neither of those are currently supported.

1. should be quite easy to add in, But i simply haven’t had the time for it atm.

2. that one i dont know how one would easily implement as it is non-linear, all degrees of freedom are connected through the integral constraint (+ having to solve 1. first)

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 5, 2024, 3:23pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/16 "2024-05-05T15:23:07Z")

</div>

Regarding 2 above, the system I want to solve is the Laplace equation, so I think determining K via the integral constraint is a linear problem. That said, without the capability to do 1 in particular, I don’t think that the problem I originally posed at the top of this thread is solvable with the current resources in fenicsx. Do you agree? This surprises me because everything is linear with inhomogeneous linear constraints. Was this solvable in the earlier fenics?

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [May 5, 2024, 5:35pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/17 "2024-05-05T17:35:07Z")

</div>

Legacy fenics had support for real spaces, Which would solver the int u dx=1 and u(x\_i)=u(x\_j) if the meshes align such that x\_i and x\_j has a dof at that location. It would silently fail of this mapping did not exist.

Dolfinx\_mpc extended this to handle u\_k = sum\_i=1,i≠k^M alpha\_i u\_i  
which covers a broad range of constraints (But not all constraints).

Note that dolfinx\_mpc is open source, and Im happy for others to contribute to it. I have alot of other projects on my plate these days, so adding new functionality to dolfinx\_mpc is not on the top of my list of priorities.

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [May 5, 2024, 6:30pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/18 "2024-05-05T18:30:56Z")

</div>

You have already been very gracious with your time, and I appreciate your help. As a new user, my focus has been to understand what the current capabilities of FEniCSx are and to see if I can use it for problems of interest to me. If I ever get to the point where I understand the code that is already up and running, I would also be interested in modifying it to meet my own needs and those of others. Right now that point seems to be very far off in the future. Thanks again for your help!

---

<div class="post-metadata">

**Author:** ![Daniel.Baker](https://avatars.discourse-cdn.com/v4/letter/d/b9e5f3/32.png) [@Daniel.Baker](https://fenicsproject.discourse.group/u/Daniel.Baker)\
**Post date:** [June 25, 2024, 4:18pm UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/19 "2024-06-25T16:18:03Z")

</div>

I am coming back to this problem now that I understand FEniCSx a little better. To implement a current condition as above, one must introduce a new variable V, which is the unknown constant voltage imposed on the boundary Gamma2. This voltage V must be determined by solving the constraint that the integral of the normal gradient of u over Gamma2 is required to be equal to 1, as noted in the condition initially posted in the problem. It seems to me that the best way to do this is to rewrite a Newton Solver where I can expand the original matrix A by an extra row and column to include the current contraint as a way to determine V.

There are some good descriptions of how to calculate the matrix A and the vector b, but I am not sure what sort of operations I would need to expand this matrix by adding an additional row and column to solve for V. Are there some examples out there that someone could point me too? So far, I have been looking at  
[https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html](https://jsdokken.com/dolfinx-tutorial/chapter4/newton-solver.html)  
as well as Section 5.3.2 of the book  
 ![image](https://global.discourse-cdn.com/free1/uploads/fenicsproject1/original/2X/d/d63a306e9031fd892e2b692d992ebe85b383ecc1.png)  
If needed, I could supply more details about how I want to do this, but perhaps what I have written above is enough for someone to point me to the right material. Thanks in advance!

---

<div class="post-metadata">

**Author:** ![raboynics](https://avatars.discourse-cdn.com/v4/letter/r/4bbf92/32.png) [@raboynics](https://fenicsproject.discourse.group/u/raboynics)\
**Post date:** [June 26, 2024, 7:46am UTC](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491/20 "2024-06-26T07:46:51Z")

</div>

Not sure if I can be of any help, but your problem looks like what I am doing in certain ways. I solve 2 coupled PDEs, one of them is Laplace equation for the voltage, where I specify a given current (the flux of voltage in some way) passing through 2 boundaries like in your case. As voltage is defined up to an arbitrary constant, I do not impose any Dirichlet b.c. for the voltage on any surface, I let the solver do its job and if it’s unable.to converge, I can deal with specifying the null space. The thing is that the function ‘‘volt’’ or electrostatics potential is determined automatically just given the input/output current.  
If you’re interested, have a look at the codes I posted in my recent post history.

[Next page](https://fenicsproject.discourse.group/t/implementing-a-potential-problem-with-a-current-condition/14491.md?page=2)
