Error when attempting to create a Dirichlet and Neumann boundary condition on cube

Errors
1

Hi, I’m having an issue when trying to apply both a Dirichlet BC at boundary 1, and a Neumann at boundary 2 (as defined in the code). I have a working example that produces an error for me, also shown below.

import fenics as fn

# Parameters
k = 10

g = fn.Expression('-10', degree=1)
T_0 = fn.Constant(20.0)

# Mesh Creation
dimx = 15
dimy = 10
dimz = 10
n = 2

mesh = fn.BoxMesh(fn.Point(0, 0, 0), fn.Point(dimx, dimy, dimz), n*dimx, n*dimy, n*dimz)

# Function Space Initialization
V = fn.FunctionSpace(mesh, "CG", 2)

# Variational Problem Definition
u = fn.TrialFunction(V)
v = fn.TestFunction(V)
f = fn.Constant(1)

# Boundary Conditions
class Left(fn.SubDomain):
    def inside(self, x, on_boundary):
        tol = 1
        return on_boundary and fn.near(x[0], 0, tol)

class Right(fn.SubDomain):
    def inside(self, x, on_boundary):
        tol = 1
        return on_boundary and fn.near(x[0], dimx, tol)

left = Left()
right = Right() 

domains = fn.MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
left.mark(domains, 1)
right.mark(domains, 2)

boundaries = fn.MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

dx = fn.Measure('dx', domain=mesh, subdomain_data=domains)
ds = fn.Measure('ds', domain=mesh, subdomain_data=boundaries)

bc0 = fn.DirichletBC(V, T_0, boundaries, 1)

# Function initialization
F = (k*fn.dot(fn.grad(u), fn.grad(v))*fn.dx(0) + k*fn.dot(fn.grad(u), fn.grad(v))*fn.dx(2)\
     - g*v*ds(2)\
     - f*v*fn.dx(0) - f*v*fn.dx(2))

a, L = fn.lhs(F), fn.rhs(F)


# Write solution to file
file = fn.File("Solutions/output_neumann_cube.pvd")

u = fn.Function(V)

fn.solve(a == L, u, bc0)

file << u

Here is the error I receive:

Traceback (most recent call last):
  File "file.py", line 70, in <module>
    fn.solve(a == L, u, bc0)
  File "file.py", line 233, in solve
    _solve_varproblem(*args, **kwargs)
  File "file.py", line 273, in _solve_varproblem
    solver.solve()
RuntimeError: 

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***     https://fenicsproject.discourse.group/
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the 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.
*** Process: 0
*** 
*** DOLFIN version: 2019.2.0.13.dev0
*** Git changeset:  ubuntu
*** -------------------------------------------------------------------------

This is my first time attempting to define boundaries using the Measure() function and it very well may be the source of the issue but I am not sure.

Thanks in advance!

2

You have no integral in f restricted to dx(1). Why is that?
If you only want to solve a problem on a subset of the grid, you have to use SubMesh or MeshView.

3

Hi Dokken, thanks for the response. I suppose that’s a good question, I am not sure how to correctly implement the integral on dx(1), would you be able to provide an example if you don’t mind?

I do think I could have been clearer with my goal, I am happy to have the problem solved across the whole surface continuously, the only thing I wanted to implement is a Neumann boundary condition only on the surface defined at 2 and a Dirichlet only on 1. Unless I’m mistaken, I only need to solve once over the entire surface, right?

4

Let’s take a few steps back.
In your code you use

and

while what you probably want to do use to just use fn.dx. What is the reason for restricting the cell integrals by integers?

5

My logic is that I may want to eventually apply different k-values in the two different domains, but for the time being, I am content with them being equal. I used the following tutorial as a general guide (https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=measure), which employs this method in 2D.

Am I mistaken in my understanding that the integer following fn.dx specifies which domain it integrates over? So I specified that both domain 0 and domain 2 as I wish for them to be both integrated over.

The only reason I had intended to designate the dx into parts is in case I want to change the k-values and that I used the boundary labels in my Dirichlet BC initialization

Edit: My goodness, I have realized I was defining dx as its own element, then I was calling fn.dx() instead of my own definition. I will update if this is indeed the solution to the issue.

Edit 2: It seems to be able to very slowly solve the problem, but I’ll have to analytically check it to see if it’s right.

6

that will partially resolve your issue, as it will then use the subdomain data supplied in:

However, I would probably use a DG-0 function to indicate different materials rather than different constants. Then you could do:

Q = FunctionSpace(mesh, "DG",0)
k = Function(Q)

# assign data to k, see and use `dx` measure.

See

A motivation for doing it this way is to avoid extra generated code, which happens often in ufl:

1 Like