Unable to solve mixed Poisson equation


39
views
0
12 days ago by
Leo  
Hi
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
11 days ago by
In your workaround, try replacing the call to solve() with
solve(A , z.vector(), b)​
This should resolve the type error.  
It solved the issue although when I want to see the results in Paraview I get this error:

Cannot read point data array "f_17-1" from PointData in piece 0.  The data array in the element may be too short.​

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

written 11 days ago by Leo  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »