Is it possible to solve such Quasi-Statics problem in Dolfin/FeniCS?

Hello,
   I think it would be constructive for you to try building the model yourself and posting codes demonstrating problems as they arise. Some of the demos will be useful for guidance:

demo_poisson.py demonstrates how to develop the necessary variational forms for the poisson equation
demo_neumann_poisson.py will help with the neumann boundary conditions
demo_submesh.py will help with working with SubDomain to work with your box inside of a box
for setting your external current density you can use a constant
J=Constant((0,0,1e9))
which should appear in your variational form. You can use a vector valued expression Expression((...,...,...)) if your current should have a spatial dependence. There are a number of other poisson related demos as well that could be of use.

Good luck,
  Nathan Borggren, PhD

Hello, Nathan!

Thank you. I tried to write some code.
Submesh and subdmain examples are very complicated, so I tried to create finite dipole. I started from demo_neumann_poisson.py example.
It does not work, it reports zero potential on big cube.

from dolfin import *

# Create mesh and define function space
mesh = UnitCube(10, 10, 10)
V = FunctionSpace(mesh, "CG", 1)
R = FunctionSpace(mesh, "R", 0)
W = V * R

# Define variational problem
(u, c) = TrialFunction(W)
(v, d) = TestFunctions(W)
# Form finite dipole
f = Expression("1e9*(x[0] == 0.4995) * (x[1] == 0.4995) * (x[2] == 0.4995) - 1e9*(x[0] == 0.5005) * (x[1] == 0.5005) * (x[2] == 0.5005)")
g = 0;
a = (inner(grad(u), grad(v)) + c*v + u*d)*dx
L = f*v*dx + g*v*ds

# Compute solution
w = Function(W)
solve(a == L, w)
(u, c) = w.split()

# Plot solution
plot(u, interactive=True)

By finite dipole do you mean that you are trying to place a charge at (0.4995,.4995,.4995) and another at (.5005,.5005,.5005)? Your form is such that a constant charge density has been placed throughout the cube, I think you are wanting localized charge sources. Perhaps try using a gaussian distribution to approximate your point sources and place them at those points.

I realized your suggestion as
f= Expression("1e9*exp(-(pow(x[0] - 0.4995, 2) + pow(x[1] - 0.4995, 2) + pow(x[2] - 0.4995, 2)) ) - 1e9*exp(-(pow(x[0] - 0.5005, 2) + pow(x[1] - 0.5005, 2) + pow(x[2] - 0.5005, 2)) )");

The surface potential is not zero, but its summetry differs from expected (it is diagnonal, not horizontal plane at z=0.5).

I would expect diagonal symmetry from your choice of locations, I think you mean:

f= Expression("1e9*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2) + pow(x[2] - 0.4995, 2)) ) - 1e9*exp(-(pow(x[0] - 0.500, 2) + pow(x[1] - 0.500, 2) + pow(x[2] - 0.5005, 2)) )")

I made small mistake, you are right. Thank you. The summetry now is OK. So we created point current point dipole source, right?
It is to possible to realize original idea with known current density in a small volume?

Yes, we created a point dipole source, but you would need to normalize it still to make the 1e9 meaningful. Try again with the SubDomains and SubMesh demos, you need to understand how to define and mark an internal boundary. Your variational form will contain terms integrated over dx(0), if 0 marks the inner volume, and terms integrated over dx(1), if marked 1 otherwise. You will define your current density much like you did with f, but just a vector expression or constant. I think if you fight with the demos, making sure you see how the forms are obtained from the differential equations, you should have no trouble solving this. This isn't really the forum to develop a problem from scratch, we need codes that break in order to fix them!

Thank you for your answers, Nathan.
I'll try to learn examples, do some trial-and-error attempts and reopen the question if I face problems.