I want to experiment with the very weak solution for Laplacian, like this simple problem
in Omega: -\Delta u = 0, in Omega,
u in L^2(\Delta, L^2(\Omega)) := {u \in L^2, \Delta u \in L^2}
on boundary: u=g, g \in H^(-1/2)(\partial\Omega)
I have made a search in the documentation and I didn’t find any reference to this simple problem.
From implementation viewpoint, we take g in P0(\partial\Omega) and look for u in P0(\Omega), satisfying
\int_\Omega (u*\Delta uh) = 0, for all uh smooth; u = g on bdry
Here, the space of u, P0, is different from the space of uh, Uh. A reasonable choice for the space of uh is Uh={uh_i, i=1, Th.nt}, with uh_i of class C^2 and with support in the union of all triangles Th_ij, where Th_ij is any triangle adjacent to the triangle Th_i.
I guess that this very weak formulation requires some proficiency with Freefem coding (especially for handling the space Uh). Can anyone help how to implement this simple algorithm in Freefem?
It seems to me that the discrete problem you want to solve is ill-posed, simply because
the space of test functions has dimension higher than that of the space in which you are looking for the solution.
What if you take a standard discrete formulation u,u_h\in P1. The question is then wether this discrete solution converges to the continuous weak solution as h\to 0.
But maybe you have already tried that as well as other standard choices of finite element spaces for some particular g which is not much regular, and found that the discrete solution does not converge?
Yes, here the space of test functions is different from the space where u lives. From this perspective the problem is ill-posed. But, this is expected from the beginning since u is discontinuous and uh are C^2.
However, we can formulate a well-posed problem from the algebraic viewpoint.
For example u in P0:
u=\sum_{i=1,Th.nt} u_i*phi_i, {phi_i, i=1, Th.nt} base of P0,
and
uh in {uh_1,…. uh_{Th.nt} a family of C^2 functions, with uh_j C^2 with compact support in
\cup_{Th_jk, Th_jk=Th_j, or a triangle adjacent to Th_j}, j=1,...,Th.nt
So, we will have a system of Th.nt equations, which a priori defines a good problem (we can impose some BC, say of Dirichlet type).
Of course, the solution - whatever it is, will depend on the choice of uh_j. But this is another problem.
Then for a square system, you need to be able to compute each coefficient of your matrix, which seems to be \int_{T_i}\Delta u_j. If you choose your u_j you should be able to compute that and build your matrix.