load cell data given for rectangular grid and use it in program

Question

I am solving the standard Poisson equation as shown in the demo, on a unit square.
My problem is that the coefficient $k(x,y)$ is given to me as piecewise constant data in the form of a $10x10$ array, where the value at (i,j) denotes the value of $k$ in the square cell located at the ith row and jth column.

How do I construct an expression $k$ from the given data that can then be used in fenics as usual to solve the Poisson equation ?

Answer 1

After some digging around and some trial and error, I found two ways to do it. Posting them both here.

nx = 3
ny = 3
lx = 1.0
ly = 1.0
mesh = RectangleMesh(0, 0, 1, 1, nx, ny, "left")
import numpy as np
k_data = np.ones(nx * ny)

# Method 1
V0 = FunctionSpace(mesh, 'DG', 0)
k  = Function(V0)

for cell_no in range(mesh.num_cells()):
  centroid = Cell(mesh, cell_no).midpoint()
  xmin, ymin = 0.0, 0.0
  deltax, deltay  = lx/nx, ly/ny
  i, j = 0, 0

  while i < nx:
      if between(centroid.x(), (xmin, xmin + deltax)):
          break
      xmin += deltax
      i +=1

  while j < ny:
      if between(centroid.y(), (ymin, ymin + deltay)):
          break
      ymin += deltay
      j += 1

  print i, j
  k.vector()[cell_no] = k_data[i+ j*nx]
## End of method 1 # k is now a function


# Method 2
k0 = MeshFunction("double", mesh, 2)
for cell in cells(mesh):
    print cell
    centroid = cell.midpoint()
    xmin, ymin = 0.0, 0.0
    delta_x, delta_y = 1.0 / nx, 1.0 / ny
    i, j = 0, 0

    while i < nx:
        if between(centroid.x(), (xmin, xmin + delta_x)):
            break
        xmin += delta_x
        i += 1

    while j < ny:
        if between(centroid.y(), (ymin, ymin + delta_y)):
            break
        ymin += delta_y
        j += 1

    k0[cell] = k_data[i + j * nx]


class MyExpression1(Expression):

    def eval_cell(self, value, x, ufc_cell):
        value[0] = k0.array()[ufc_cell.index]

k = MyExpression1()
## End of method 2 # k is now an Expression