This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

how to make my form work correctly, one order equations

0 votes

i have a problem as follow
\nabla p = f
where p is a scalar and f is 2\times 1 vector
the boundary condition of p is p_bc
I want to know whether my Fenics form is right
my form: inner(grad(p), V)*dx - inner(f,V) = 0

THE SYSTEM TOLD ME THAT:
*** -------------------------------------------------------------------------
*** Error: Unable to set given (local) rows to identity matrix.
*** Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where: This error was encountered inside PETScMatrix.cpp.
*** Process: 0

COULD YOU TELL ME HOW TO MAKE MY FORM CORRECT.

asked Sep 4, 2016 by Hamilton FEniCS Novice (500 points)

F2 = inner(grad(p),v)dx - inner(project_R,v)dx
a2 = lhs(F2)
L2 = rhs(F2)

solve(a2 == L2, p1, p_bc)

Hi, please provide the minimal complete code the shows your problem.

from dolfin import *

solve problem \nabla p = f with dirichlet boundary conditions ......

Create mesh and define function space

mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
Mh = FunctionSpace(mesh, "Lagrange", 1)

Define Dirichlet boundary (x = 0 or x = 1)

def boundary(x):
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS

Define boundary condition

u0 = Constant(0.0)

u0 = Expression(["0.0","0.0"])
bc = DirichletBC(V, u0, boundary)

Define variational problem

the function spaces are right ?????

f_grad = Expression(["2picos(2pix[0])sin(2pix[1])","2pisin(2pix[0])cos(2pix[1])"])

Define variational problem

p = TestFunctions(Mh)
v = TrialFunctions(V)
f = Constant((0, 0))
a = inner(grad(p),v)dx
L = inner(f,v)
dx

Compute solution

u = Function(Mh)
solve(a == L, u, bc)

Save solution in VTK format

file = File("poisson.pvd")

file << u

Plot solution

plot(u, interactive=True)

The way you write it, the problem is ill-posed. For illustration, consider first the problem without the boundary conditions

from dolfin import *

# solve problem \nabla p = f

mesh = UnitSquareMesh(32, 32)
V = VectorFunctionSpace(mesh, "Lagrange", 1)
P = FunctionSpace(mesh, "Lagrange", 1)

p = TrialFunction(P)
v = TestFunction(V)

a = inner(grad(p), v)*dx

A = assemble(a)
print 'A is %d x %d' % (A.size(0), A.size(1))

You are getting a rectangular matrix. Consider instead defining $p$ as the solution to the minimization problem of $|\nabla q - f|_{L^2}$ for $q\in P$.

...