Solving a Poisson type equation

variational formulation
1

Hello,

I am trying to solve the pde \nabla^2 \phi = -\nabla^2 \psi for the potential \phi in a 3D domain, which is insulated except at one point where it is earthed. The function \psi is known on the 3D mesh and consists of a solid region of 1’s with the remainder being 0. My variational formulation is as follows:

def boundary_fixed_point(x,on_boundary):
    tol=1E-10
    return (abs(x[0]-3.0) < tol) and (abs(x[1]-0.0)<tol) and (abs(x[2])<tol)

V = FunctionSpace(mesh,'CG',1)
phi = TrialFunction(V)
psi = TrialFunction(V)
v = TestFunction(V)
psi = Function(V)
# I read in the pointwise definition of psi from a file. psi is defined at the  nodes of the mesh 
V_g = VectorFunctionSpace(mesh,"CG",1)
grad_psi = project(grad(psi),V_g)
dbc = DirichletBC(V,Constant(0.0),boundary_fixed_point,
    method='pointwise')
a = inner(grad(phi),grad(v))*dx
L = inner(Constant(-1.0)*grad_psi,grad(v))*dx

phi = Function(V)

solve(a == L,phi,dbc)

While this code code does run, it produces incorrect answers. The answers produced do not agree with another code using the same mesh and input conditions, nor the physics of the problem in really simple problems. I assume that I must have either the variational formulation incorrect or there is some problem with the calculation of grad_psi.

If someone could provide some insight that would be great. I can expand on my explanation of the problem or the code if necessary.

Thank you,

Peter.

2

Hi,
I know its a while since this post was initially created, but I figured better late than never. You are missing several RHS terms stemming from integration by parts.
I’ve tried to recreate your problem, using a DG0 function to create two domains.

from dolfin import *
mesh = UnitCubeMesh(17,17,17)
def boundary_fixed_point(x,on_boundary):
    tol=1E-10
    return (abs(x[0]-0.5) < tol) and (abs(x[1])<tol) and (abs(x[2])<tol)

V = FunctionSpace(mesh, "CG", 1)
x = SpatialCoordinate(mesh)
psi = 0.1*cos(pi*x[0])*cos(pi*x[1])*cos(pi*x[2])
Q = FunctionSpace(mesh, "DG",0)
delta = Function(Q)
coords = Q.tabulate_dof_coordinates()
for i, coord in enumerate(coords):
    if coord[0] <= 0.5 and 0.5 <= coord[1] and coord[2] <= 0.5:
        delta.vector()[i] = 1
u, v = TrialFunction(V), TestFunction(V)
n = FacetNormal(mesh)
a = inner(grad(u), grad(v))*dx
l = -delta*inner(grad(psi), grad(v))*dx\
    +jump(delta)*inner(jump(grad(psi)), n("+"))*jump(v)*dS\
    +delta*inner(grad(psi), n)*v*ds

uh = Function(V)
bc = DirichletBC(V,Constant(0.0),boundary_fixed_point,
                 method='pointwise')
solve(a==l, uh, bc)
XDMFFile("uh.xdmf").write(uh)
XDMFFile("delta.xdmf").write(delta)