# 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. – Wolfgang Bangerth Feb 5 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 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. – Abdullah Ali Sivas Feb 5 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 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. – Abdullah Ali Sivas Feb 7 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.$$
• 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$. – Wolfgang Bangerth Feb 6 at 17:47