### Unable to solve mixed Poisson equation

39

views

0

Hi

I am trying to solve a mixed Poisson equation. Here is the code:

I am trying to solve a mixed Poisson equation. Here is the code:

```
from dolfin import *
import numpy as np
from meshfile import mesh2
mesh = mesh2
class BL_Bottom(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.00488 and x[0] <= 0.00613 and (x[1] >= 0.00488 ) and (x[1] < 0.005 + DOLFIN_EPS ) else False
class BL_Left(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.00488 and x[0] < 0.005 + DOLFIN_EPS and (x[1] >= 0.00488 ) and (x[1] <= 0.01013 ) else False
class BL_Top(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.00488 and x[0] <= 0.00613 and (x[1] > 0.01 - DOLFIN_EPS ) and (x[1] <= 0.01013 ) else False
class BL_Right(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] > 0.006 - DOLFIN_EPS and x[0] <= 0.00613 and (x[1] >= 0.00488 ) and (x[1] <= 0.01013 ) else False
def boundary(x):
return x[0] == 0.00488 or x[0] == 0.00613 or x[1] == 0.00488 or x[1] == 0.01013
domains = CellFunction("size_t", mesh)
domains.set_all(0)
B_L_Bottom = BL_Bottom()
B_L_Bottom.mark(domains, 2)
B_L_Left = BL_Left()
B_L_Left.mark(domains, 2)
B_L_Top = BL_Top()
B_L_Top.mark(domains, 2)
B_L_Right = BL_Right()
B_L_Right.mark(domains, 2)
#lable domain
dx = Measure('dx', domain=mesh, subdomain_data=domains)
#Defining element for scalar variables (e.g. concentrations and voltage)
Element1 = FiniteElement("CG", mesh.ufl_cell(), 1)
#Defining element for vector variable (e.g. displacement)
Element2 = VectorElement("CG", mesh.ufl_cell(), 2)
# Defining the mixed function space
W_elem = MixedElement([Element1, Element1 ])
W = FunctionSpace(mesh, W_elem)
#c_ref_plus, c_ref_minus = split(z)
c_ref_plus, c_ref_minus = TrialFunctions(W)
#Defining the test functions
(v_3,v_4) = TestFunctions(W)
# Define Dirichlet boundary conditions for the in BL for c_ref_plus
bc_1 = DirichletBC(W.sub(0),Constant(0.0), boundary)
bc_2 = DirichletBC(W.sub(1),Constant(0.0), boundary)
#Collecting the boundary conditions
bcs = [ bc_1,bc_2]
#Summing up variational forms
a = dot(grad(c_ref_plus), grad(v_3)) * dx(2) + dot(grad(c_ref_minus), grad(v_4)) * dx(2)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)
L =f*v_3*dx(2) + f*v_4*dx(2) + g*v_3*ds(2) + g*v_4*ds(2)
f = File("P.pvd")
z = Function(W)
#solver
solve(a == L, z, bcs)
(c_ref_plus,c_ref_minus) = z.split(True)
f << c_ref_minus
```

I get this error :```
Error: Unable to set given (local) rows to identity matrix.
Reason: some diagonal elements not preallocated (try assembler option keep_diagonal).
```

I also tried this for the solver:

```
A = assemble(a, keep_diagonal=True)
b = assemble(L)
# Apply boundary conditions
[bc.apply(A, b) for bc in bcs]
solve(A , z, b)
(c_ref_plus,c_ref_minus) = z.split(True)
```

But it did not work! It gives me this error:

`TypeError: in method 'la_solve', argument 2 of type 'dolfin::GenericVector &'`

Any idea how to solve this issue?

Thanks!

Community: FEniCS Project

### 1 Answer

1

In your workaround, try replacing the call to solve() with

`solve(A , z.vector(), b)`

This should resolve the type error.
Please login to add an answer/comment or follow this question.

Do you know what is wrong with that and how I can visualize the results?