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  
It looks to me that this example may help you

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  

1 Answer


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)
Please login to add an answer/comment or follow this question.

Similar posts:
Search »