Poisson Problem: Unable to set given (local) rows to identity matrix.

Question

Hi,

When I run the following Poisson Problem, I get an 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.

UFL:

element = FiniteElement("Lagrange", triangle, 1)

u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
g = Coefficient(element)
D1 = Coefficient(element)
D2 = Coefficient(element)

a = inner(D1*grad(u), grad(v))*dx(1) + inner(D2*grad(u), grad(v))*dx(2)
L = f*v*dx

main.cpp

class left_side : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
      return(x[1] <= 0.5);
  }
};

class right_side : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  {
      return(x[1] >= 0.5);
  }
};

class DirichletBoundary : public SubDomain
{
  bool inside(const Array<double>& x, bool on_boundary) const
  { return on_boundary; }
};

int main()
{
  UnitSquareMesh mesh(32, 32, "crossed");
  // Create mesh functions over the cell facets
  MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim());
  subdomains = 0;
// MeshFunction<std::size_t> subdomains(mesh, mesh.topology().dim(), 0);
// CellFunction<std::size_t>subdomains(mesh, (dolfin::MPI::rank(mesh->mpi_comm()) == 0));

  left_side Left;
  Left.mark(subdomains, 1);
  right_side Right;
  Right.mark(subdomains, 2);

  Poisson::FunctionSpace V(mesh);
  DG0::FunctionSpace Vdg(mesh);
  Constant f(1.0);
  Constant D1(1.0);
  Constant D2(10.0);
  Constant u0(0.0);
  DirichletBoundary boundary;
  DirichletBC bc(V, u0, boundary);

  Poisson::BilinearForm a(V, V);
  Poisson::LinearForm L(V);
  a.D1 = D1; a.D2 = D2;
  L.f = f;

  Function u(V);
  solve(a == L, u, bc);
  return 0;
}

Answer 1

Hi,
I think that you need to redefine the measure dx as (in python, feel free to look for the c++ equivalent):

      dx = Measure("dx")(subdomain_data=subdomains)

Otherwise, dx is not aware of subdomains. The default is to have an unique subdomain with value 0. Thus, when you assemble the linear form in your code, dx(1) and dx(2) loop over empty subdomains, resulting in a zero matrix without preallocated elements on the diagonal.

Best regards,
Francesco

Comments

Thanks for your reply Francesco. That's the issue. I have been trying to do that in C++ but have not been able to. I couldn't find any example in undocumented/documented demos or anywhere else either.

---

I have not tested it, but thanks to
https://fenicsproject.org/qa/10310/subdomains-in-c
I think it should be enough for you to modify your code as follows:

   ....
   a.D1 = D1; a.D2 = D2;
   a.dx = subdomains
   ....

---

Thank you so much! It works perfectly. Such a simple solution and I was trying all complex things.