Solving for scalar function from its gradient

Question

I would like to know how to quickly find a scalar function $a$ given the values of $\partial a / \partial x$ and $\partial a / \partial y$, with the scalar function having a known value at a point.

Two solutions I'm working on are the following:
Solving an inverse problem using derivative() within fenics or using dolfin-adjoint
If I want just one $a$ value at a point I could integrate from the known point to the required point on a straight line

Has anyone needed to do this calculation before? Are there any functions in fenics that can do this?

Answer 1

I figured it out, its a linear problem, I can do it as follows:

V = VectorFunctionSpace(mesh, 'P', 2)
Q = FunctionSpace(mesh, 'P', 1)

PPoint  = 'near(x[0], 0.0) && near(x[1], 0.0)'
bcp = DirichletBC(Q, Constant(0.0), PPoint, 'pointwise')

p  = TrialFunction(Q)
v = TestFunction(V)
F5 = dot(as_vector([px_,py_]) - grad(p), v)*dx
a  = lhs(F5)
L  = rhs(F5)

p = Function(Q)
solve(a == L, p, bcs=bcp)

But I get a dolfin error:

*** Error:   Unable to apply boundary condition.
*** Reason:  Dofmap ownership range (0,1089) does not match matrix row range (0,8450).
*** Where:   This error was encountered inside BoundaryCondition.cpp.
*** Process: 0

The domain is a 0 to 1 square. What is the reason for this error?

Comments

Hi, the error is due to the assembled matrix being rectangular.

A = assemble(a)
print A.size(0), A.size(1)

---

Thank you for the answer Miro. Is the assembled matrix rectangular because of the function spaces I'm using? Shouldn't the assembly matrix always be square?

---

I got the gradient inverse by solving a Poisson equation as described in Efficient Computation of the Inverse Gradient on Irregular Domains. The code is shown below:

PPoint  = 'near(x[0], 0.0) && near(x[1], 0.0)'
bc      = DirichletBC(Q, Constant(0.0), PPoint, 'pointwise')

# Define variational problem
p = TrialFunction(Q)
q = TestFunction(Q)

f = div(as_vector([px_,py_]))
n = FacetNormal(mesh)
g = dot(n,as_vector([px_,py_]))

a = -inner(grad(p), grad(q))*dx
L =  f*q*dx + g*q*ds

# Compute solution
p = Function(Q)
solve(a == L, p, bc)

I'm still not sure why my original solution doesn't work

---

The matrix you assemble from a bilinear form has as many rows as is the dimensionality of the function space from which the test functions come and as many columns as the dimensionality of the space of trial functions. Think of the block off-diagonal block inner(div(v), q)*dx in the Stokes problem for example.

---

Thank you for the reply Miro. I changed the order of the V functions space to 1 to try to match the number of rows but I still get the same error:

*** Error:   Unable to apply boundary condition.
*** Reason:  Dofmap ownership range (0,1089) does not match matrix row range (0,2178).
*** Where:   This error was encountered inside BoundaryCondition.cpp.
*** Process: 0

The order of the px_ and py_ functions is 1 as well. The error message says that the error is from BoundaryCondition.cpp, could it have something to do with how I'm writing the boundary condition? Could the problem be my use of the as_vector() function?

The original pde I'm trying to solve for in 2D is:

$$p - \nabla u = 0$$

Where u is the p function in the code and p is composed of the two scalar functions px_ and py_.

---

The issue is that you are still assembling a rectangular matrix.