Assembled matrix is zero

1

I’m trying to solve Poisson equations on the layered cylindrical mesh with Fenics 2019.1.0.

subdomains
subdomains905×602 45.2 KB

For each subdomain there is a conductivity tensor basis_conds defined. For some values of basis_conds, however, the assembled matrix contains only zeros (A.nnz() = 0). Here is the code that I’m using:

mesh = femg.fwd_fem.read_mesh_xdmf(os.path.join(data_dir, '4shell_cylinder_mesh.xdmf'))
subdomains = femg.fwd_fem.read_mesh_function_xdmf(mesh, os.path.join(data_dir, '4shell_cylinder_subdomains.xdmf'), 3)

basis_conds = {8: 1. * np.eye(3),
               9: 2. * np.eye(3),
               10: 3. * np.eye(3),
               11: 4. * np.eye(3)}
V = df.FunctionSpace(mesh, 'CG', 1)
u = df.TrialFunction(V)
v = df.TestFunction(V)

# Describe bilinear form as a som of integrals over subdomains
dx = df.Measure("dx")(domain=mesh, subdomain_data=subdomains)
form_list = []
for label in basis_conds.keys():
    cond_tensor = df.as_matrix(basis_conds[label])
    form_list.append(df.inner(df.dot(cond_tensor, df.grad(u)),
                              df.grad(v)) * dx(label))
a = sum(form_list)
A = df.assemble(a)

If I use multipliers 1, 2, 3, 4 for basis_conds, A is well assembled. However, if for example I use multipliers 8, 9, 10, 11, matrix A contains only zeros. If I define these multipliers as random positive floats between 0 and 1, some times A is well assembled, some times not.

Where lies the problem and what should I change to make the algorithm robust to variation in conductivity tensors?
Thank you in advance.

2

are you talking about the dict key or the value multiplied by eye(3)?

3

About the value multiplied by eye(3)

4

I guess the numpy eye matrix has default type as an integer. Change it by doing np.eye(3,dtype=np.float64)

5

But I multiply it by a float so the resulting array is float. And why it works for some multipliers, and not for others?

6

Only if you change the dtype. As you haven’t supplied a minimal working example code that reproduces the behavior, there isn’t much I can do to help you debug this. Please create a minimal working code that reproduces this (that is a code that I can run on my computer that produces the error).

7

Can I upload mesh files somewhere that you can run the code?

P.S. I checked, all arrays are in the dict are floats

8

Please make a code Where you do not have to use a custom mesh. Is there a particular value that doesn’t work? If so, integrate it over the whole domain of for instance a unit cube and check if your matrix is non Zero.

9

I couldn’t reproduce the error using simple cube mesh. I did however just change
(1) df.dot(cond_tensor, df.grad(u)) by
(2) cond_tensor * df.grad(u)
and it looks like it is now working.
To conclude, I don’t know what is the problem but the fact is that for the same data and input using (1) results in empty assembled matrix while using (2) works well.