Variable coefficient in LHS

Question

Hi,
I'm pretty new to Fenics, so the solution to my question may be obvious, but I'm still struggling to crack on it. Basically, I have a Coefficient in the LHS of my variational equation that varies along the Grid nodes:

# Equation
a = ((k**2 / n**2) * inner(E, v) - inner(nabla_grad(E), nabla_grad(v))) * dx 

That's the 'k' coefficient and it has different values along the grid. I have all these values stored in a matrix of the same order as the Mesh grid (100x100 grid and a 100x100 matrix with all the k values).

Is there a simple way to make the solver take account of such spatial variable coefficient values?

Thanks in advance,
Patricio,

Answer 1

You might define k before:
Example:

def k(x,y):
   return x+y
a = ...k(x[0],x[1])...

Does this work?

Comments

Well, in my particular scenario "k" isn't an expression but a coefficient that takes different values across the space - it's a refraction index particularly.

This way, for example a region of the grid may belong to air or to a concrete wall and have a different refraction index than their adjacent node....That's why I store refraction index information in a Matrix.

---

Hi, so you are looking for a method to represent your table data as a CG1 function, correct? Also, is this for a structured mesh and what is the logic behind the table? I mean is table[i, j] the value of index at some point [ih_x, jh_y]?

---

The table has no logic, it depends on spatial information gotten from a .png file. The .png file describes refraction indexes in a certain room for instance. So the matrix will contain refraction indexes for externa, intermediate walls, doors, etc.
That means that one matrix element may be 2.55+0.5i (concrete refraction index) and the adjacent one just 1+0i (air refraction index).

Sorry for not being clear enough on that in my first note.

Answer 2

Hi

The variable coefficient does not usually pose any additional difficulties at for a steady simulation if you put the values of k in a Function. So it seems to me like your problem is really how to get k from you matrix (structured mesh?) to a proper Function? This is a bit difficult since Fenics uses an unstructured mesh, and I suggest you search for answers including vertex_to_dof_map, dof_to_vertex_map. To get a structured matrix k onto an unstructured UnitSquareMesh, though, you can try something similar to the following:

from dolfin import *
import numpy

N = 10
mesh = UnitSquareMesh(N-1, N-1)
V = FunctionSpace(mesh, 'CG', 1)

# To get indices in the structured mesh reverse-engineer using x = i/N, y = j/N
x = interpolate(Expression("x[0]*{}".format(N-1)), V)
y = interpolate(Expression("x[1]*{}".format(N-1)), V)

# Create some NxN matrix of k values:
ks = numpy.outer(numpy.arange(N), numpy.arange(N))

# Put ks in k
k = Function(V)
v2d = vertex_to_dof_map(V)
for d in v2d:
    ii = int(x.vector()[d])
    jj = int(y.vector()[d])
    k.vector()[d] = ks[ii, jj]

plot(k)

Comments

Mikael,
That sounds like a good option. To provide more detail about the matrix, I basically read a .png file (once read it becomes a matrix structure) and based on the pixel colors I map them to refraction indexes - that result is also a matrix. You're right, my problem was to map Mesh unstructured indexes to matrix elements. I'll apply your suggestions and circle back with an official outcome,

Thanks!

---

Mikael,
That sounds like a good option. To provide more detail about the matrix, I basically read a .png file (once read it becomes a matrix structure) and based on the pixel colors I map them to refraction indexes - that result is also a matrix. You're right, my problem was to map Mesh unstructured indexes to matrix elements. I'll apply your suggestions and circle back with an official outcome,

Thanks!

****UPDATE***** - this option has worked fine, I've just had to modify this line this way inside the loop:

k.vector()[d] = ks[jj, ii]

I'm getting message not related to this after I run my script:

/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py:1448: ComplexWarning: Casting complex values to real discards the imaginary part
vec_values *= values

And the solution field I get is completely real (no Im components). k is complex-valued as well. Is there any way to get around this?

Thanks!

---

Also, when I switch from a UnitSquareMesh to a RectangleMesh I start to get indexing issues ......Is there some different in the indexing between these 2 mesh types? If so, what would be an appropiate correction?

/usr/lib/python2.7/dist-packages/dolfin/cpp/la.py:1448: ComplexWarning: Casting complex values to real discards the imaginary part
vec_values *= values
Traceback (most recent call last):
File "", line 4, in
IndexError: index 360 is out of bounds for axis 1 with size 360