### Construct least-squares system

215
views
0
4 months ago by
I have a 2D mesh and some data points (x_i, y_i), z_i in the mesh. (x_i, y_i are not necessarily node points.) I would now like to solve Laplace's equation on the mesh with natural boundary conditions such that the resulting piecewise function fulfills u(x_i, y_i) = z_i as closely as possible. This can be done by formulating a linear least-squares problem à la
$- \Delta u = 0 \quad\text{in } \Omega,\\ n \cdot \nabla u = 0 \quad\text{in } \partial\Omega,\\ u(x_i, y_i) = z_i \quad\forall i.$
I could use https://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.sparse.linalg.lsqr.html for solving the above problem, but I'll have to give it the system matrix.

Any hint on how to construct it using FEniCS? Or is there perhaps even a least-squares solver in FEniCS itself?
Community: FEniCS Project
Getting the system matrix is easy. But I dont know how to construct your constraint equation. It would be of the form B*u=z where entries of B are the shape functions evaluated at your data points.
written 4 months ago by Praveen C

https://fenicsproject.org/olddocs/dolfin/1.4.0/python/demo/documented/neumann-poisson/python/documentation.html

The only modifications that I would do are
f=Constant(0)
g=Constant(0)​
written 4 months ago by Ruben Gonzalez
Hi Nico, your question is a bit unclear regarding what you mean by a LS solver. A priori, your solutions are in Sobolev spaces, thus pointwise operations are meaningless. Also, you can get the matrix of the assembled system without the points z_i, but it still would have the problem of having a non trivial kernel without even looking at your point z_i. If it was a LS problem where you want to impose the PDE, then your problem would require a Lagrange multiplier somewhere. In general, it doesn't seem quite well posed as it standes, and maybe we need more information.

Best regards!
written 4 months ago by Nicolás Barnafi
Thanks for the replies!

> Hi Nico, your question is a bit unclear regarding what you mean by a LS solver

I'm looking at the equation system after discretization of course; sorry if that was unclear. A better way of writing it down is perhaps as the equivalent norm-minimization problem.

> A priori, your solutions are in Sobolev spaces, thus pointwise operations are meaningless.

Well, it's a bit like Dirichlet boundary conditions; those are meaningless too. I guess the common thing to do is to introduce traces, but no one ever talks about that anymore. I silently ignored that here, too; thanks for pointing out the inconsistency though.
written 4 months ago by Nico Schlömer

3
4 months ago by
The Laplacian part of the equation can be assembled as usual; not quite as trivial is the  u(x_i, y_i) part. As answered in [this question](https://www.allanswered.com/post/lkbkm/evaluating-all-trial-functions-at-a-handful-of-points/), one can use evaluate_basis_all à la
bbt = BoundingBoxTree()
bbt.build(mesh)
dofmap = self.V.dofmap()
el = self.V.element()
for xy in centers:
cell_id = bbt.compute_first_entity_collision(Point(*xy))
cell = Cell(mesh, cell_id)
coordinate_dofs = cell.get_vertex_coordinates()

v = numpy.zeros(1, dtype=float)
el.evaluate_basis_all(
v, xy, coordinate_dofs, cell_id
)

print(v[0])​

to construct the evaluation matrix.

The matrix can then be converted into a scipy csr matrix via

# [...] create PETScMatrix()

row_ptr, col_indices, data = L.mat().getValuesCSR()
size = L.mat().getSize()
A = sparse.csr_matrix((data, col_indices, row_ptr), shape=size)