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