# L2 projection with bounds

In some problems we are currently working on, we are working with discontinuous functions that are defined on a finite element mesh and are established using Lagrangian particles. To obtain them on a finite element basis, we calculate the L2 projection. This works quite well and gives us a piecewise continuous representation of these functions.

One issue that we face is that we would like to force our L2 projection to be within given bounds (without loss of generality, say [0,1]).

Is there a way to force an L2 projection to be within a certain set of bounds?

Note: We solve for the L2 projection by solving the following finite element problem on $$\Omega$$. Where $$\phi$$ are Lagrangian polynomials (e.g. Q1).

$$\int_\Omega \phi_i \sum_j \phi_j \epsilon_j = \int_\Omega \phi_i f$$

• Are you willing to restrict yourself to (bi-/tri-)linear functions on each cell? It is simple to achieve what you want in that case, but practically impossible if you are interesting in going to quadratic functions. Feb 5, 2021 at 17:44
• Let's assume that bi/tri linear functions are sufficient for what we need :). I can foresee that at higher order it is not the same cup of tea in terms of difficulty.
– BlaB
Feb 5, 2021 at 19:14
• I feel like you are not asking about your main problem, but rather you are asking about your approach to solve your main problem. I think the solution to your problem is flux limiters or slope limiters which prevent overshoots and undershoots when you use a DG-like method with explicit time steppers to solve a hyperbolic PDE. Feb 5, 2021 at 21:29
• No this is exactly my main problem. I need to define a field using lagrangian particle and a L2 projection and I would like to maintain it between bounds. This not related to fluxes or DG methods...
– BlaB
Feb 5, 2021 at 21:32
• I see. Flux limiters, I believe, are still applicable. Especially if your finite element solution is piecewise polynomial (i.e., if $u\in\text{DG}$). But Wolfgang Bangerth's solution is valid and probably, much easier as an approach. Feb 7, 2021 at 17:38

Let's stick with the case of either piecewise constants or (bi-/tri-)linear functions onto which to project first, and I'm going to talk about "linear" in the following even though on quadrilateral/hexahedral meshes, that's not quite correct. In any case, you are trying to solve the following problem if you have no constraints on the values of your interpolant: $$\min_{u \in X} \frac 12 \sum_i \left(u(x_i)-v_i\right)^2$$ Here, $$X$$ is the set of all linear functions, and $$v_i$$ are the values associated with your particles, which are located at positions $$x_i$$. This is a standard least squares problem and, because you can write $$u(x) = \sum_j U_j \varphi_j(x)$$ is equivalent to finding coefficients $$U_i$$ so that you minimize $$\min_{U \in {\mathbb R}^n} \frac 12 \sum_i \left(\sum_j U_j \varphi_j (x_i)-v_i\right)^2.$$ In other words, it's a quadratic minimization problem, nothing difficult here.

If, now, you want to restrict the set you're searching in to $$X_{[0,1]} = \{u \in X: 0\le u(x) \le 1 \; \forall x\}$$ then the first question you'd have to ask is whether $$X_{[0,1]}$$ is a convex set. The answer to this is yes, so there is hope that this is a well-behaved optimization problem, and in fact it is relatively easy to solve because we can characterize these functions well: Since (bi-/tri-)linear functions attain their minima and maxima at the boundary of the the cell, you have that $$X_{[0,1]} = \{u \in X: u(x)=\sum_j U_j \varphi_j(x), 0\le U_j \le 1 \}$$ So it remains a quadratic optimization problem in $$n$$ variables, but now with $$n$$ bound constraints on the variables. This is, again, not very difficult to solve. One could use a simple active-set method, for example, and would probably be done with an effort two or three or so times the cost of just the unconstrained optimization problem. Just to write it out, the problem you're now trying to solve is of the form $$\min_{U \in {\mathbb R}^n} \frac 12 U^T M U - U^T F \\ \qquad \text{such that } 0 \le U_j \le 1, \qquad 1\le j\le n.$$

Things do become substantially more expensive and, in practice, mostly intractable if you used (bi-/tri-)quadratic or even higher order shape functions. That's because for quadratic functions, the minimum and maximum of a function is not attained at the node points: You can have a function with all positive expansion coefficients and still the function can be negative in some places. So, in general, $$\max_x u(x) \neq \max_j U_j$$ and the same for the minimum. In other words, if you wanted to write $$\min_{u \in X_{[0,1]}} \frac 12 \sum_i \left(u(x_i)-v_i\right)^2$$ in terms of the expansion coefficients, this is not the same as $$\min_{U \in {\mathbb R}^n} \frac 12 U^T M U - U^T F \\ \qquad \text{such that } 0 \le U_j \le 1, \qquad 1\le j\le n.$$

To the best of my knowledge, there is no easy way to characterize which coefficients lead to a function whose values are between 0 and 1. If you did want to do things like this, you'd have to have a basis with shape functions that are non-negative (e.g., the Bernstein basis), and that would solve the issue with the lower bound: non-negative coefficients lead to a non-negative function. But I doubt that you could find a basis that could satisfy both the constraints you want to place on your projection.

As a consequence, I don't think there is an easy way to achieve enforcing bounds on quadratic functions. From all I know, people only do this for linear functions in practice.

• Thanks for the amazing answer. If I understand correctly. If you were to solve the same problem but using quadratic basis function, the consequence would be that your nodal values would be between 0 and 1, but since these function allow for local minima within cells, you would be left with no guarantee that you don't exceed 1 or go below 0 within your elements. Am I correct? If that is so, I can live with that :)!
– BlaB
Feb 5, 2021 at 21:26
• Precisely. You can see that in 1d. Look at the quadratic shape function associated with one of the two vertices. So your nodal values are $(1,0,0)$, but the function is then $\phi_1(x)=2(1-x)(\tfrac 12 - x)$ which has a minimum of $-1/8$. Feb 6, 2021 at 17:47