Boundary conditions in mixed formulation

Question

Hi,

I have a trouble with the definition of boundary conditions in mixed formulation :

from dolfin import *
mesh = UnitSquareMesh(10, 10)
facet_f = FacetFunction('size_t', mesh, 0)
CompiledSubDomain('near(x[0], 0)').mark(facet_f, 2)

ds = Measure("ds", domain=mesh, subdomain_data=facet_f)

P1 = VectorFunctionSpace(mesh, "Lagrange", 1) 
P2 = FunctionSpace(mesh, "Lagrange",1)
V = P1 * P2

mu, Kpar = 1., 10.

bc1 = DirichletBC(V.sub(0).sub(1), Constant(0.0), facet_f, 2)
bcs = [bc1]

calT = Expression(("t2"),t2=0.)

w = Function(V)
(u,p) = w.split();

I = Identity(2)
FF = I + grad(u)
C = FF.T*FF
Ic = tr(C)
JJ  = det(FF)

psi = (mu/2.*(Ic-2) - p*(JJ-1))*dx
W = psi + calT*u[0]*my_ds(2) 
v = TestFunction(V)
du = TrialFunction(V)

F = derivative(W,w,v)
J = derivative(F,w,du)

problem = NonlinearVariationalProblem(F,w, bcs, J=J)
solver  = NonlinearVariationalSolver(problem)
solver.parameters.nonlinear_solver = "snes"
solver.parameters.snes_solver.linear_solver = "umfpack"
solver.solve()

I obtain the following error :

*** 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.

This error disappears when I suppress the boundary condition in my code....

Do someone have an idea about it ?

many thanks in advance !

Claire

Comments

Hi, I only suggest a workaround for the issue as I don't see a way how to pass the keep_diagonal flag to the assembler through the current API of the Nonlinear* classes. The idea is to introduce a "do nothing" modification to the problem. The following seems to work

W = psi + calT*u[0]*ds(2) + Constant(0)*inner(w, w)*dx

---

I am not sure to understand why, but if I replace :

(u,p) = w.split()

by :

(u,p) = split(w)

it works !

FEniCS syntax can be tricky !