Using data-matrix as source term in possion equation

Hi everyone,

I am trying to solve a Poisson equation with a right-hand side that depends on data that I measured and is stored in a csv-file/matrix. Basically an expression like:
‘\Laplace (n) = U(x,y)’,
where U is a scalar field and n the function I am trying to get. That leads me to the following weak formulation to solve:
‘\int \nabla(m) \nabla(n) dA = \int U m dA’,
where m is my testfunction. (Assuming Neumann-BCs for now)
I also have the x and y data stored in separate matrices and I use them to generate a corresponding mesh.
I managed to import the data and generate a MeshFunction out of it. To use it within the variational problem, I tried to adopt the tensor-weighted poisson example to generate an expression in C++. However, when trying to solve it, I get the following error:

 UFLException: Can only integrate scalar expressions. 
 The integrand is a tensorexpression with value shape (1,) and free indices with labels ().

So I guess something is wrong with the way i am trying to use the expression.
My code looks something like that (simply using some artificial data here):

import numpy as np
import csv
import dolfin as df
import matplotlib.pyplot as plt

#Data import
U = np.ones([122,208])
#Mesh generation
[nx, ny] = U.shape
lx =  16.4485;
ly =  28.1388;
x_min = -8.4508;
y_min = -6.4118;
x_max = 7.9977;
y_max = 21.727;
rect_mesh = df.RectangleMesh(df.Point(x_min,y_min), df.Point(x_max,y_max), nx, ny, "left")

#Rearange the imported data to a Mesh-fitting format:
U0 = df.MeshFunction("double", rect_mesh, 2)
for cell in df.cells(rect_mesh):
    #print(cell)
    centroid = cell.midpoint()
    xmin, ymin = x_min, y_min
    delta_x, delta_y = lx / nx, ly / ny
    i, j = 0, 0

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

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

    U0[cell] = U[-i, -j]

    ##################################


# Code for C++ evaluation of displacement U:
displacement_code = """

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;

#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>

class Displacement : public dolfin::Expression
{
public:

  // Create expression with 1 components
  Displacement() : dolfin::Expression(1) {}

  // Function for evaluating expression on each cell
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
  {
    const uint cell_index = cell.index;
    values[0] = (*U0)[cell_index];
  }

  // The data stored in mesh functions
  std::shared_ptr<dolfin::MeshFunction<double>> U0;
};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<Displacement, std::shared_ptr<Displacement>, dolfin::Expression>
    (m, "Displacement",py::module_local())
    .def(py::init<>())
    .def_readwrite("U0", &Displacement::U0);
}

"""

U0_expression = df.CompiledExpression(df.compile_cpp_code(displacement_code).Displacement(), U0=U0, degree=0)
##################################
    
# Compute solution
P1 = df.FiniteElement("Lagrange", rect_mesh.ufl_cell(), 1)
R = df.FiniteElement("Real", rect_mesh.ufl_cell(), 0)
W = df.FunctionSpace(rect_mesh, P1 * R)
    
# Define variational problem
(n, c) = df.TrialFunction(W)
(m, d) = df.TestFunctions(W)
rhs = (U0_expression*m)*df.dx;
lhs = -(df.inner(df.grad(n), df.grad(m)))*df.dx
#L =  f*v*df.dx - g*v*df.ds

W = df.Function(W)
df.solve(lhs == rhs, W)

# Save solution in VTK format
file = df.File("poisson.pvd")
file << u

# Plot solution
plt.plot(u)
plt.show()

I am using FEnics in anaconda / spyder. My Dolfin version is 2019.1.0.
If anyone has a solution to solve this issue or another approach to use external data in such way, I’d be happy to read them! Thanks in advance
moritz

The simplest solution seems to be just commenting-out the constructor in the C++ code:

  // Create expression with 1 components
  //Displacement() : dolfin::Expression(1) {}

Thanks, seems to be working!