### Construct least-squares system

215

views

0

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?

\[

- \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

### 1 Answer

3

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)
```

Please login to add an answer/comment or follow this question.

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

The only modifications that I would do are

Best regards!

> 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.