How to use a numpy matrix to define the conductivity

Hi everyone,

I’m trying to use Fenics to solve the diffusion equation: grad(*kgrad(u)) = 0**

In my situation, k is a scalar field defined from a given numpy matrix M as following:
k = Expression(‘M[int(x[0]*nx/l),int(x[1]*ny/h)’, M=M, nx=nx, ny=ny, l=l, h=h)

This is not working as M is not recognized by C++ I guess.

Can anyone help me to solve this issue?

Kind regards

Welcome to the community!

Your code shouldn’t compile. For one, there’s the missing closing square bracket in M[... and you should need to supply the either the degree or the element argument to Expression.

There’s a recent thread where @kamensky posted the most recent solution on how to pass a numpy array to an Expression:

You may be able to circumvent the hassle by defining k as a Function and manually setting the values at the DOFs as described here:
http://home.simula.no/~hpl/homepage/fenics-tutorial/release-1.0-nonabla/webm/materials.html

1 Like

Hi Klunkean,

Thank you!

I found these two lines of code which enable you to assign a numpy array M to a Dolfi vector in a simple manner. Each element of M corresponds to the physical value of a given scalar field (permeability for instance ) at each node.

k = Function(V)
k.vector()[:] = M

M is a numpy array of shape (n, ) with n denotes the number of elements.

Regards

This is of course assuming the ordering of your matrix elements corresponds to the dof ordering in k.
And from the information you supplied I assumed your matrix had 2 columns and as many rows as you have elements in a 1D domain. Then your method actually shouldn’t work straight away. But glad it works for you.

Also consider using DG0 elements for the function space of k if you want cell-wise constant values in k.

Good day Klunkean,

I’m glad to hear from you.

This is perfectly working for me. I actually get the coordinates of the elements using the command:

V.tabulate_dof_coordinates()

Then I use these coordinates to convert my initial 2D matrix into a 1D vector corresponding to the same ordering of K.

I will have a look on DG0. I feel it’s doing the same thing I’m doing, but in a simpler fashion :smiley:

Have a good day!