Non-diffusive Variational Problems with Distributional and Weak Gradient Constraints

In this paper, we consider non-diffusive variational problems with mixed boundary conditions and (distributional and weak) gradient constraints. The upper bound in the constraint is either a function or a Borel measure, leading to the state space being a Sobolev one or the space of functions of bounded variation. We address existence and uniqueness of the model under low regularity assumptions, and rigorously identify its Fenchel pre-dual problem. The latter in some cases is posed on a non-standard space of Borel measures with square integrable divergences. We also establish existence and uniqueness of solutions to this pre-dual problem under some assumptions. We conclude the paper by introducing a mixed finite-element method to solve the primal-dual system. The numerical examples confirm our theoretical findings.

1. Introduction Some Bibliography 1.1.Organization of the paper.2. Notation and Preliminaries 2.1.The gradient constraint 3. Existence Theory for (P) 3.1.The case when α is a nonnegative Borel measure 3.2.The case when α is an integrable function 4. Existence Theory for the pre-dual problem (P * ) 4.1.The case when p is a function and α is either a function or a measure 4.2.The case when α is a function and p is a measure 5. Duality relation between (P) and (P * ) 5.1.The case when α is a function 5.2.The case when α is a measure 6.A Finite Element Method with Applications 6.1.Finite Element Discretization 6.2.Numerical Examples References

Introduction
We begin by considering an evolutionary problem whose semi-discretization (in time) gives rise to the class of stationary problems of interest in this paper.Suppose that f : (0, T ) × Ω → R together with u 0 : Ω → R are given, where Ω ⊂ R d is a bounded domain with a Lipschitz boundary.Furthermore, let α be a given nonnegative function (possibly only integrable), or a nonnegative Borel measure in Ω. Suppose that u : (0, T ) × Ω → R, such that u(0) = u 0 , is a solution to the following problem where the set K is given by K := Ũ (0, T ) ∩ {w : w(t) ∈ K almost everywhere}. (1. 2) The choice of Ũ (0, T ) and K in (1.2) hinges on the type of the boundary conditions and the regularity of α.We assume that the boundary ∂Ω is partitioned into a Dirichlet boundary part Γ D and a non-Dirichlet boundary part Γ N , both composed of a finite number of connected parts, such that Notice that on Γ N , we do not necessarily prescribe Neumann boundary conditions, as it is later clarified.However, a conservation law of material is in place in the case Γ D = ∅; specifically, it can be inferred from (1.1) that ´Ω(u(T ) − u 0 ) dx = ´T 0 ´Ω f dx dt given that v = u ± 1 are admissible test functions as we see next.The restriction of u to the Γ D part of the boundary is assumed to be zero, and no restrictions are assumed on Γ N .The set K is convex and it arises by a nonlinear law with a bound on the first order derivative terms.In the most general form K is given by with 1 ≤ p ≤ +∞.We briefly discuss the two possible scenarios that we consider: (I) If α is a nonnegative measurable function, then U Γ D (Ω) is a Sobolev-type space and G = ∇ is the weak gradient, so that |∇v| p is the p -norm of the weak gradient of v. Hence, |∇v| p ≤ α in (1.3) is considered in the almost everywhere (a.e.) in Ω sense.(II) If α is a nonnegative Borel measure, then U Γ D (Ω) is a subset of functions of bounded variation BV(Ω).In this case, G = D is the distributional gradient, and |Dv| p the total variation measure of Dv associated to the p -norm, and the constraint |Dv| p ≤ α is understood in the measure sense.Both instances, (I) and (II) are related, in fact (I) may be considered as a special case of (II).Furthermore, letting α ∈ M + (Ω) in case (II), where M + (Ω) denotes the set of nonnegative Borel measures, enables us to handle the delicate case α ∈ L 1 (Ω) + in (I).Next we shall provide a brief description of modeling capabilities of (I) and (II) in the context of a particular application.
A possible motivation for the above class of problems is based on the study of accumulation of granular heterogeneous material on possibly discontinuous structures.This approach was pioneered by Prigozhin [28,30,31] in the case of homogeneous materials and a continuous support structure.In this vein, f : (0, T ) × Ω → R represents the (density) rate of a granular material being deposited on a supporting structure u 0 : Ω → R.Moreover, ´T 0 ´Ω f dx dt is the total amount of material deposited on Ω over the time interval [0, T ].In case that α > 0 is a real number, this corresponds to the classical case of a granular cohesionless material where homogeneous piles are generated.If α : Ω → R is not constant zero, the value of α at a point determines the angle of repose of the < l a t e x i t s h a 1 _ b a s e 6 4 = " P n x A u 8 9 3 Z 2 c W 6 q a G 3 S X u m l g F 4 Z k = " > A A A B 7 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 W P R i 8 c K 9 g P a U C b b T b t 0 s w m 7 G 6 G E / g g v H h T x 6 u / x 5 r 9 x 0 + a g r Q 8 G H u / N M D M v S A T X x n W / n d L a + s b m V n m 7 s r O 7 t 3 9 Q P T x q 6 z h V l L V o L G L V D V A z w S V r G W 4 E 6 y a K Y R Q I 1 g k m d 7 n f e W J K 8 1 g + m m n C / A h H k o e c o r F S p 4 8 i G W N l U K 2 5 d X c O s k q 8 g t S g Q H N Q / e o P Y 5 p G T B o q U O u e 5 y b G z 1 A Z T g W b V f q p Z g n S C Y 5 Y z 1 K J E d N + N j 9 3 R s 6 s M i R h r G x J Q + b q 7 4 k M I 6 2 n U W A 7 I z R j v e z l 4 n 9 e L z X h j Z 9 x m a S G S b p Y F K a C m J j k v 5 M h V 4 w a M b U E q e L 2 V k L H q J A a m 1 A e g r f 8 8 i p p X 9 S 9 q 7 r 7 c F l r 3 B Z x l O E E T u E c P L i G B t x D E 1 p A Y Q L P 8 A p v T u K 8 O O / O x 6 K 1 5 B Q z x / A H z u c P w 0 u P M A = = < / l a t e x i t > ↵ < l a t e x i t s h a 1 _ b a s e 6 4 = " d j R h g 0 G 2 L U A F c h 7 e e P 0 V D k n m i D o = " > A A A B 6 H i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E 0 W P R i 8 c W 7 A e 0 o W y 2 k 3 b t Z h N 2 N 0 I J / Q V e P C j i 1 Z / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / + g f H j U 0 n G q G D Z Z L G L V C a h G w S U 2 D T c C O 4 l C G g U C 2 8 H 4 b u a 3 n 1 B p H s s H M 0 n Q j + h Q 8 p A z a q z U C P v l i l t 1 5 y C r x M t J B X L U + + W v 3 i B m a Y T S M E G 1 7 n p u Y v y M K s O Z w G m p l 2 p M K B v T I X Y t l T R C 7 W f z Q 6 f k z C o D E s b K l j R k r v 6 e y G i k 9 S Q K b G d E z U g v e z P x P 6 + b m v D G z 7 h M U o O S L R a F q S A m J r O v y Y A r Z E Z M L K F M c X s r Y S O q K D M 2 m 5 I N w V t + e Z W 0 L q r e V d V t X F Z q t 3 k c R T < l a t e x i t s h a 1 _ b a s e 6 4 = " N G q 5 x J Y g 8 Q t 6 Z P 8 U e 8 A T e y 4 q n y U r e 2 d 4 q 7 p b 3 9 g 8 O j 8 v F J R 0 e J I r R N I h 6 p X o A 1 5 U z S t m G G 0 1 6 s K B Y B p 9 1 g e p / 5 3 S e q N I t k y 8 x i 6 g s 8 l i x k B J t M S q q t y 2 G 5 Accumulation of two kinds (magenta and blue) of granular materials on discontinuous surface.(LEFT) Depiction of f (t, x) = f 1 χ (x 0 ,x 2 ) (x) + f 2 χ (x 2 ,1) (x), the accumulation of both materials, and dα = α The value of the initial supporting structure u 0 and the the final distribution u(T ).material at that point, i.e., the steepness of a cone generated from a point source of material.This is the case for heterogenous sandpiles [7] and also a restricted case of the quasi-variational sandpile model; see [3,4,27,5].In a more general setting, where α is a measure, using the approach in this paper, it is possible to generate discontinuous structures such as cliffs by preserving discontinuities in the initial supporting structure u 0 and/or of f .Such an approach has not yet been considered by the literature to the best of our knowledge.
A description of the qualitative behavior of Problem (1.1) is displayed in Figure 1.We assume two materials with different angles of repose α 1 and α 2 with α 1 > α 2 are poured on the discontinuous structure u 0 (x) := χ (x 1 ,x 2 ) (x) for x ∈ Ω := (0, 1) and 0 < x 1 < x 2 < 1.The intensity of the material being deposited is given by f (t, x) = f 1 χ (x 0 ,x 2 ) (x) + f 2 χ (x 2 ,1) (x) for some points x 0 , and x 2 , and some f 1 , f 2 > 0, i.e., the first and second materials are poured with density rates f 1 and f 2 , respectively, during the entire time interval (0, T ).We further assume that a sharp edge can form at x 2 with maximum height of 1, and in addition discontinuities of maximum size 1 can be preserved at the locations of the discontinuities of u 0 .Finally, the the gradient constraint α is then given by dα = α 1 χ (0,x 2 ) (x) dx + α 2 χ (x 2 ,1) (x) dx + 3 i=1 δ(x − x i ), and the material is assumed to escape freely at the boundary points of Ω.On the right side of Figure 1, we see the comparison between u 0 and u(T ), the solution at time T > 0; on the left we see the depiction of f , α, and the accumulation regions of both materials.
The study of solutions to (1.1) usually makes use of the semi-discretization (in time) of the problem via an implicit Euler method.In particular, we approximate the partial time derivative ∂ t u by (u n − u n−1 )/k for some time-step k > 0. The class of problems of interest in this paper is then given by Minimize (min) Notice that (P) can be seen as the time-discrete version of (1.1) where the solution u to (P) is equal to u n when f corresponds to ´nk (n−1)k f (τ ) dτ + ku n−1 with u n−1 given.Closely related to the problem above, we consider the following class of problems We prove that (P * ) is the Fenchel pre-dual of problem (P), i.e., the Fenchel dual [14] of (P * ) under certain conditions is (P).Several choices for V Γ N (Ω) and J are explored which are directly related to the nature of α.In all cases considered, V Γ N (Ω) contains d-dimensional vector fields with divergences in L 2 (Ω).In particular, we consider (i) If α is a nonnegative measurable function (additional assumptions are later explained but continuity is enough to guarantee what follows), then we explore two options for J: In the first case Borel measures, so that the second functional denotes the integral of α with respect to the total variation measure of p induced by the q -norm.The two functionals are closely related, and the first can be seen as a restriction of the second one to measures that are absolutely continuous with respect to the Lebesgue measure.(ii) If α is a nonnegative Borel measure, then V Γ D (Ω) is contained in the space of maps that are α measurable, with J given by A few words are in order concerning (P) and (P * ).Although the objective functional in (P) is smooth and amenable, the constraint set K makes the entire problem highly nonlinear and nonsmooth.The latter also holds for (P * ) given the nature of the functional J.The development of solution algorithms for both problems is a rather delicate issue that requires appropriate regularization methods that can handle the nonsmothness in an asymptotic fashion.
The paper focuses on functional analytic properties of (P) and (P * ) together with duality relationship properties.Additionally, we develop a mixed finite type method to solve the optimality conditions corresponding to (P) and (P * ).Some Bibliography.The structure of Problems (P * ) and (P) and their inherent difficulties are analogous to the ones that appear in the context of plasticity; see [25,35] and references therein.In particular, the first class of applications for diffusive variational problems with gradient constraints is the elasto-plastic torsion problem.Such a problem has been thoroughly analyzed by Brézis, Caffarelli, Evans, Friedman, Gerhardt, and others; see [16,15,8,18,17,10,11,9]. Further, a complete account of the literature can be found in [33].A significant amount of the aforementioned works focuses on regularity of solutions, the free boundary, and the equivalence of the gradient constrained problem to a double obstacle one.
The modeling of the evolution of the magnetic field in critical-state models of type-II superconductors also leads to a problem like (1.1) with the addition of a diffusive operator and a state-dependent constraint in some cases; see [32,2,29,30,20,23,21].See [34] for a study of evolutionary variational problems with non-constant gradient constraints, and [26] for a complete account of evolutionary problems with derivative bounds.
Analogous problems are found in mathematical imaging involving total variation regularization [24,19,6] and more specifically in the weighted total variation version [22].There, in contrast to the work here, the L ∞ -norm on the gradient is replaced by the L 1 -norm, leading to a pre-dual problem with a pointwise bound in its state variable.
1.1.Organization of the paper.Preliminaries are provided, and some notations are made explicit in Section 2; elementary results about the generalized gradient constraint are given in Section 2.1.In Section 3, we prove existence and uniqueness of the solution to problem (P) for the cases when α is either a nonnegative Lebesgue measurable function or a nonnegative Borel measure.Existence of solutions to problem (P * ) is addressed in Section 4, while for the case when p is a function we require d = 1, when p is a measure the dimension restriction is dropped.The relation between problems (P) and (P * ) are considered in Section 5, where a rigorous Fenchel duality result establishes a link between these two problems.In particular, in Section 5.1, we address the case where α is a function and the variable p is either a function or a measure.The case when α is a measure and an extension of the duality result of the previous section is given in Section 5.2.Finally in Section 6, we introduce a mixed finite element method to solve the underlying problems and present a range of numerical tests.

Notation and Preliminaries
The purpose of this section is to introduce notation involving spaces, and convergence notions that are used throughout the paper; in particular, we address the well-known notions of Sobolev spaces and the space of functions of bounded variation.We refer the reader to Attouch et al. [12] that we follow closely for this introduction together with the book of Adams and Fournier [1].
For a Banach space X, we denote its corresponding norm as • X .For an element F in the topological dual X of X, the duality pairing of F and an arbitrary element x ∈ X is written as F, x X ,X .Throughout the paper, all Banach spaces are assumed to be real vector spaces.
The inner product on the Lebesgue space L 2 (Ω) of (equivalence classes of) functions that are square integrable on Ω is denoted as (•, •), so that (f, g) := ´Ω f (x)g(x) dx for f, g ∈ L 2 (Ω) where dx refers to integration with respect to the Lebesgue measure.
The Sobolev space of functions in L r (Ω) for 1 ≤ r < +∞ with weak gradients in L r (Ω) d is denoted by W 1,r (Ω), and it is endowed with the norm where ∇v denotes the weak gradient of v.In the case r = 2, we use the notation H 1 (Ω) := W 1,2 (Ω).Given that Ω is assumed Lipschitz, restriction of a function v ∈ W 1,r (Ω) to the boundary ∂Ω is welldefined via the continuous trace map γ 0 : W 1,r (Ω) → L r (∂Ω).Furthermore, the closed subspace of functions in W 1,r (Ω) that are zero on Γ D is denoted by . The space of real-valued Borel measures M(Ω) is endowed with the norm µ M(Ω) := |µ|(Ω), where |µ| is defined for an arbitrary open set O as Note that µ, z M(Ω),C 0 (Ω) = ´Ω z dµ, and that |µ| defines a Borel measure in M + (Ω), the subset of non-negative elements of M(Ω), i.e., σ ∈ M + (Ω) if σ(B) ≥ 0 for every Borel set B ⊂ Ω.
We denote by BV(Ω), the space of functions v in L 1 (Ω), for which the total variation semi-norm is finite and where q is the Hölder conjugate of p, i.e., 1/p + 1/q = 1; see [12,Section 10.1].The space BV(Ω) is a Banach space endowed with the norm The operator D represents the distributional gradient, and for a v ∈ BV(Ω), Dv is a R d -valued Borel measure.We use |Dv| p to denote the total variation measure (associated to the p -norm) of Dv, and the total mass |Dv| p (Ω) is by definition Furthermore, the Lebesgue decomposition result applied to Dv implies that there exist measures D a v and D s v such that Dv = D a v + D s v, with D a v and D s v respectively being absolutely continuous and singular with respect to the ddimensional Lebesgue measure.
We define now the notions of weak and intermediate convergence of sequences in BV(Ω) which provide different topologies on the space BV(Ω).The former is obtained by a subsequence of a bounded sequence in BV(Ω).Moreover, the latter is sufficient to preserve boundary conditions in the sense of the trace as stated in Theorem 2.3 below.Definition 2.1 (Weak convergence for BV(Ω)).Let {u n } be a sequence in BV(Ω) and u * ∈ BV(Ω).We say that u n converges to u * weakly, denoted as u n u * in BV(Ω), if The definition 2.1 is understood in light of the following fact: If {u n } is a bounded sequence in BV(Ω), there exists u * ∈ BV(Ω) such that along a subsequence u n u * in BV(Ω).The latter follows since the embedding BV(Ω) → L 1 (Ω) is compact (see Attouch et al. [12,Theorem 10.1.4.]) for Lipschitz domains, and since a bounded sequence of measures admits a weakly convergent subsequence.
We shall use the direct method of calculus of variations to establish existence of solutions to problems in BV(Ω) and with Dirichlet homogeneous boundary conditions on Γ D .The space of interest is BV Γ D (Ω) defined as where γ 0 is a trace operator; see [12, section 10.2].Notice that we use the same notation for the trace operator in Sobolev spaces W 1,p (Ω).There is a fundamental issue with the trace in BV(Ω) and the application of the direct method as we show next with an standard example adapted from [12].
Consider a bounded sequence {u n } in BV Γ D (Ω).Then, we can extract a subsequence (not relabeled) of {u n } such that u n u * in BV(Ω) .The problem is that in general it is not possible to say that u * ∈ BV Γ D (Ω): Let Ω = (0, 1) with Γ D = {0}, and consider {v n } defined as Then, v n ∈ BV Γ D (Ω), and The underlying reason is that the trace operator in BV(Ω) is not continuous with respect to weak convergence, but it is with respect to the intermediate convergence subsequently defined.We further notice that |Dv n |(0, 1) = 1 and |Dv * |(0, 1) = 0, this discrepancy is central to the issue we are considering.Definition 2.2 (Intermediate convergence).Let {u n } be a sequence in BV(Ω) and u * ∈ BV(Ω).We say that u n converges to u * in the sense of intermediate convergence if and The name intermediate convergence arises since it describes a stronger topology than the one of weak convergence, but not as strong as the norm one.The importance of the intermediate convergence can be seen in the following result which holds in our case since Ω ⊂ R d is a Lipschitz bounded domain.We refer to Attouch et al. [12,Theorem 10.2.2] for its proof.2.1.The gradient constraint.A few words are in order concerning the gradient constraint given in the set K defined in (1.3).Although in the case when G = ∇ the situation is somewhat standard, if G = D, the distributional gradient for BV functions, require several non-trivial explanations.In the cases where α is a Borel measure and v ∈ BV(Ω), the inequality Since the inequality in (2.5) holds for every w n by initial assumption, it also holds in the limit for w.Furthermore, (2.5) immediately implies (2.2), so the result is proven.

Existence Theory for (P)
In this section, we discuss the existence and uniqueness of solution to the problem (P).We start with the case when α is a measure, and the case when α is a function follows as a special one.In particular, existence of solutions is studied in the function spaces Both of these spaces share the same difficulty: Bounded sequences do not necessarily admit convergent (in some sense) subsequences that preserve the zero boundary condition on Γ D in the limit.The main purpose of this section is to overcome this obstacle.
3.1.The case when α is a nonnegative Borel measure.We consider in this section that α ∈ M + (Ω) and hence the state space is given by We start by proving the following lemma which gives sequential precompactness of some classes of bounded sets in BV Γ D (Ω).These bounded sets are subsets of K which in this case is defined as The above results particularly means that for a sequence {v n } in K that is bounded in BV(Ω), there exists a subsequence that converges to some u * ∈ BV(Ω) in the sense of intermediate convergence.Further, u * ∈ BV Γ D (Ω) and also u * ∈ K.A direct consequence of the above lemma is the following result.Theorem 3.2.If α ∈ M + (Ω), then there exists a unique solution to (P) in BV Γ D (Ω).
Proof.Consider an infimizing sequence {u n } for (P).It follows that {u n } is bounded in L 2 (Ω) and hence Lemma 3.1 is applicable.That is, there is a subsequence of {u n } (not relabelled) such that u n u * in L 2 (Ω), and u n → u * in the sense of the intermediate convergence for BV(Ω), and further u * ∈ K. Finally, by exploiting the weakly lower semicontinuity property of the objective functional in (P), we have that u * ∈ K is a minimizer.
Next we discuss the case when α is a function.

3.2.
The case when α is an integrable function.In this section, we let α : Ω → R be a nonnegative and integrable function, leading to Γ D (Ω).This case can be interpreted (to some extent) as a special case of the one in the previous subsection under the assumption that α is a measure absolutely continuous with respect to the Lebesgue measure.However, we proceed in a slightly different fashion by considering α as a function and the state space contained in W 1,1 (Ω); this provides further insight on bounded sequences in K and in Sobolev spaces.In this case, we have K given by Next we state a version of Lemma 3.1 adapted to the current setting which can be used to prove existence of solutions to (P).Lemma 3.3.Let α ∈ L 1 (Ω) + and M > 0, then every sequence {v n } in the set and for some v * ∈ K * , which is also the weak limit in W 1,1 Γ D (Ω) of the same subsequence.The above can be seen as a consequence of equi-integrability of the set K. Recall that a family of functions F ⊂ L 1 (Ω) is equi-integrable provided that for every > 0, there exists a δ > 0 such that for every set A ⊂ Ω with |A| < δ we have that ´A |u| dx < for all u ∈ F. Further, the Dunford-Pettis theorem states that if {u n } is a bounded sequence in L 1 (Ω) and is equi-integrable, then u n u along a subsequence for some u ∈ L 1 (Ω).Hence, since K is bounded in W 1,1 (Ω), and the gradients are equi-integrable, it is simple to infer strong convergence in L 1 (Ω) together with weak convergence of the gradients in L 1 (Ω).The improvement of the latter convergence is done again via Prokhorov's result as in the proof of Lemma 3.1 leading to an equivalent of the intermediate convergence in BV(Ω).The trace preservation follows directly from the same proof.Further note that the convergence determined does not imply strong convergence in W 1,1 (Ω) since this space is not uniformly convex.Another formulation of the above lemma is that bounded sets with equi-integrable gradients are compact in W 1,1 Γ D (Ω) when endowed with the metric With the use of Lemma 3.3 and following the same argument as before for Theorem 3.2, we have Theorem 3.4.If α ∈ L 1 (Ω) + , then there exists a unique solution to (P) in W 1,1 Γ D (Ω).

4.
Existence Theory for the pre-dual problem (P * ) The focus of this section is on existence and uniqueness of solutions of problem (P * ) under different functional analytic settings.In particular, we focus on two cases where p is either (i) a function or (ii) a Borel measure.In the first case, we let α be either a function or a measure; here, existence results are limited to d = 1.On the other hand, in the second case we establish an existence and uniqueness result for p with arbitrary d ∈ N, for a specific class of α's (to be specified later).Furthermore, this second case requires a nonstandard space of vector measures with divergences in L 2 (Ω).Remarkably, a version of the integration-by-parts formula still holds in this general setting; such a construct is rather recent [37].We start with the case when p is a function.
4.1.The case when p is a function and α is either a function or a measure.We begin this section by considering that α ∈ L 1 (Ω) + and J is defined as Moreover, we define We assume that if d = 1 and Γ N = ∅ then α is not identically zero, and if d > 1 then α > 0 a.e. in Ω.Thus, the space V Γ N (Ω) is defined by where which follows from the fact that J(p) + 1 2 ´Ω |p (x)| 2 dx is an equivalent norm (to the usual one) on H 1 Γ N (Ω).The latter is due to J(p) = ´Ω α(x)|p(x)| dx being a seminorm in H 1 Γ N (Ω) and norm on the constants, i.e. for a ∈ R, J(a) = |a|α(Ω) = 0 iff a = 0; see [36,Chapter 1.4].We can now establish existence of a solution to problem (P * ).Theorem 4.1.Let d = 1, α ∈ L 1 (Ω) + , and if Γ N = ∅ then suppose that α is not identically zero.Consider J as defined in (4.1) on V Γ N (Ω) as in (4.2).Then, there exists a unique solution to (P * ).
Proof.The proof is based on the direct method.Let J : V Γ N (Ω) → R be the objective function in (P * ), that is, and let , and there exists a weakly convergent (not relabeled) subsequence ) we have existence of a subsequence (not relabeled) p n → p in C(Ω).Finally, weak lower semicontinuity of J (p) yields that p ∈ V Γ N (Ω) is a solution to (P * ).The strict convexity of the objective functional provides uniqueness to the solution.
An analogous approach can be considered when α is a non negative Borel measure (and not identically zero), that is, when α ∈ M + (Ω).In particular, we set and we construct the space V Γ N (Ω) in the same way as in (4.2), but with the norm • α,2 defined as p α,2 := ˆΩ |p| q dα + div p L 2 (Ω) , and assuming that if d = 1 and Γ N = ∅ then α is not identically zero, and if d > 1 then α(B) > 0 if |B| > 0 and B ⊂ Ω is a Borel set.The existence result of Theorem 4.1 follows mutatis mutandis: Since 1 2 ´Ω |p (x)| 2 dx + ´Ω |p| dα is again a norm in H 1 Γ N (Ω), see [36,Chapter 1.4], the exact argument is applicable in this case.We can now focus on the case when p is a measure which provides a general setting for the problem of interest in terms of existence, uniqueness, and duality results.

4.2.
The case when α is a function and p is a measure.We focus now on problem (P * ) when J is defined as and p is a Borel measure.Notice that the above functional can be seen as a generalization of the functional in (4.1).The latter can be obtained by letting p be absolutely continuous with respect to the Lebesgue measure.
The functional analytic setting in this section, requires p to be a measure with divergence of p in L 2 (Ω), and α to be measurable with respect to |p| q .We start with a proper definition of such spaces and their properties.We disregard the possible "boundary conditions" for the variable p, so that Γ N = ∅ and we define V Γ N (Ω) as follows: where M(Ω) and we define div p := h.The space W (Ω) is a Banach space when endowed with the norm where q ∈ [1, +∞] and Note that above •, • is the duality pairing between M(Ω) d and C c (Ω) d , and hence Similarly to the definition of |w| q (Ω), we can define |w| q (A) for any open set A, and subsequently for an arbitrary Borel set A. Hence, |w| q induces a nonnegative measure (the total variation measure of w); in addition |w| q (Ω) = ´Ω d|w| q .Note that the space W (Ω) contains regular maps, clearly if p ∈ C 1 c (Ω) d then p ∈ W (Ω), in this case "d|p| q = |p| q dx" where dx is the Lebesgue measure.A note on the space W (Ω) is in order.Although one may be inclined to think that vector fields whose divergences are in L 2 (Ω) would always have better regularity than just the measure type, this is not true.We consider an example developed by Šilhavý [37] to show otherwise.Let u ∈ BV(Ω) with Ω ⊂ R 2 , and define p = (Du) ⊥ with (a 1 , a 2 ) ⊥ = (a 2 , −a 1 ) with Du the distributional (measure valued) gradient of u; it follows that div p = 0.This can be seen as follows: C ∞ (Ω) is dense (in the sense of the intermediate convergence) in BV(Ω), this means in particular that lim ´Ω ∇ϕ • p n dx = ´Ω ∇ϕ • dp for such a smooth sequence defined as p n = (Du n ) ⊥ with u n ∈ C ∞ (Ω).Since also ´Ω ∇ϕ • p n dx = 0, the result follows by taking the limit and from (4.6).
Following Šilhavý [37], we have a form of integration-by-parts formula together with a trace result.We denote by Lip B (Λ) the space of Lipschitz maps h : Λ → R for Λ ⊂ R k and endow it with the norm where Lip(h) is the Lipschitz constant of h on Λ.It follows that for each p ∈ W there exists a linear functional N p : Lip B (∂Ω) → R such that for all v ∈ Lip B (Ω) we have Further, N p is bounded in the following sense for some C > 0, and all p ∈ W and all g ∈ Lip B (∂Ω).Provided that p and v have enough differential regularity, we observe as expected.Thus, (4.8) is an extension of the usual integration-by-parts formula.
We are now ready to state and prove the existence and uniqueness result for problem (P * ) under the setting above.Theorem 4.2.Let α ∈ C(Ω) be such that α(x) > 0 for all x ∈ Ω, and consider J defined by (4.4) on V Γ N (Ω) = W (Ω) as given in (4.5).Then, problem (P * ) admits a unique solution.
Proof.Note first that J is well-defined given that α is measurable with respect to all Borel measures.Consider an infimizing sequence {p n }.Since min x∈Ω α(x) > 0, then {p n } is bounded in V Γ N (Ω).Hence, we can extract a subsequence (not relabelled) such that p n p * in M(Ω) d for some p * ∈ M(Ω) d and div Since the map p → |p| q is weakly lower semicontinuous, αp n αp * in M(Ω) d , and |q| q = α|p| q for q = αp, we have that p * is a minimizer by a weakly lower semicontinuity argument.Uniqueness follows from the strict convexity of the objective functional.
At this point, one would be tempted to extend the result to the case where Γ N = ∅, for example, by defining While the space above is well-defined, it is not clear if the weak limits of sequences in the space also belong to it.In fact, if for each v ∈ Lip B Γ D (Ω).However, the weak limit along a subsequence argument is not enough to pass to the limit in the left hand side given that ∇v is not necessarily of compact support.This remains an open problem.

Duality relation between (P) and (P * )
In this section, we discuss the dual problem corresponding to (P * ).We start with the case when α is a Lebesgue measurable function and further subdivide it into two subsections.In Section 5.1 we discuss the case when the pre-dual variable p is a function and in the following Section 5.1.2we assume that the variable p is a measure.Next in Section 5.2, we consider the case where α is a measure and the pre-dual variable p is a function.In general, we prove that Problem (P) is the Fenchel dual of Problem (P * ).
In order to keep the discussion self-contained, we introduce the following notation and terminology.For an extended real valued function ψ : X → R ∪ {∞} over a Banach space X, by ψ * we denote its convex conjugate, which is defined by (e.g.see [14, p. 16]) Provided that the operator div : V → L 2 (Ω) is defined for a Banach space V , and it is bounded, its adjoint (div) * : L 2 (Ω) → V * is well-defined and is given by (div) * v, p V * ,V = (v, div p) for all v ∈ L 2 (Ω) and all p ∈ V .

5.1.
The case when α is a function.We first consider the case where α is a non negative Lebesgue measurable function and we accordingly set in (P * ) for the cases when p is a function or a measure, respectively.For each of the choices of J above, we will also establish the strong duality to (P).We assume throughout this section (and for the sake of simplicity) that α ∈ C(Ω), and α(x) > 0, for all x ∈ Ω as discussed in Section 1, together with and hence, e. in Ω}.Note that in Section 3 we proved the existence and uniqueness of solution to (P).
We compute the dual problem to (P * ) and show that it is given by problem (P).Defining the problem (P * ) can be written as inf for div : , where the space V Γ N (Ω) is chosen based on whether p is a function or a measure.By [14, p. 61], the Fenchel dual of (P * ) with respect to the perturbation function is given by inf where the convex conjugates J * : (V Γ N (Ω)) * → R ∪ {∞}, F * : L 2 (Ω) → R ∪ {∞} of J and F are defined according to (5.1), see also [14, p. 17] for more details.
5.1.1.Duality when p is a function.Now we show that the problem (P) is the dual to problem (P * ).
In this section, we assume that V Γ N (Ω) is given by (4.2), and that We start by proving the following result: We break the proof of the above theorem into Lemmas 5.3 and 5.4, which we state after the following observation.
Lemma 5.3.Let u ∈ L 2 (Ω) with J * (div * u) = 0. Then the following hold true: Then, by using definition of a function of bounded variation, see [12, Definition 10.1.1],we have that the supremum on the right hand side of the above inequality is +∞ if u / ∈ BV(Ω) and hence, u ∈ BV(Ω) if J * (div * u) < +∞.(ii) As u ∈ BV(Ω), we have that Du ∈ M(Ω) d and the inequality |Du| p ≤ α is understood in the sense of (2.4).Hence, if for an arbitrary open set O ⊂ Ω, then the required condition |Du| p ≤ α immediately follows.By the assumption J * (div * u) = 0, and using integration by parts, we observe that where the last inequality follows using the definition of ´O |Du| p and (5.7).(iii) By (i) and (ii), it holds that ˆS |Du| p ≤ ˆS α(x) dx, ( for every Borel set S (see (2.3)), and especially for every Borel set of Lebesgue measure zero, it follows that |Du| p vanishes on every set of measure zero, and hence Du is absolutely continuous w.r.t. the d-dimensional Lebesgue measure, and therefore Du = ∇u, i.e., the distributional gradient is a weak gradient.Thus, u ∈ W 1,1 (Ω).(iv) To obtain the boundary conditions on u, we will show that if J * (div * u) = 0, then γ 0 (u) = 0.
Since u ∈ BV(Ω), then using [12, Theorem 10.2.2] we have that γ 0 (u) ∈ L 1 (∂Ω) and Subsequently for all p ∈ V Γ N (Ω) ∩ C 1 (Ω), we arrive at To get (5.9), for a p Now for ε > 0, by inner regularity [13, pp.95, proposition 15.1], there exist closed subsets Then, by Urysohn's lemma there exists Next, using the two inequalities above in conjunction with and Now since q ∈ C 1 (Ω) and ε > 0 have been chosen arbitrarily, it follows that This immediately leads to the required result, γ 0 (u) = 0 a.e. on Γ D , and the proof is complete.
Proof.Since u ∈ K, therefore by the definition of K, it holds that u ∈ W 1,1 Γ D (Ω) and |∇u| p ≤ α a.e. in Ω. Next, using the definition of the convex conjugate J * of J, we obtain that (5.10) Next, by using the density of C 1 (Ω) d ∩ V Γ N (Ω) in V Γ N (Ω), from (5.10), we obtain that Thus, since J * (div * u) is non-negative (we can set p ≡ 0 in the definition of J * ), it follows that J * (div * u) = 0 and the proof is complete.
Next we compute the conjugate function of the function F .
Proposition 5.5.The conjugate function of F defined in (5.2) is given by Proof.The proof is an immediate consequence of the definition of F * .Recalling the definition of F * and rearranging the terms, we obtain that The result then follows from elementary calculus.where u denotes the solution of (P), and ∂J(q) denotes the subdifferential of J at a point q.
Proof.As a corollary to Theorems 5.1 and 5.5, it immediately follows that the dual of problem (P * ) which is given in (5.4) is identical to problem (P).Using that J and F are convex and continuous proper functions and bounded from below, equality (5.12) and the extremality relation (5.13)  φ(p, 0) and the application of [14, (4.20) in p. 61] yields that φ is convex, l.s.c., proper, and bounded from below, given that the same holds true for J and F .Thus, it follows from [14, p. 49] that φ * * = φ and that (P * ) is identical to its bidual problem φ * * (p, 0), i.e., to the dual problem to (P), with respect to the perturbation function φ * .
Note that though we assume α ∈ C(Ω), results within this section hold for α ∈ L 1 (Ω).However recall that the existence result for this case (c.f.Section 4.1) only stands in the case d = 1.5.1.2.Duality when p is a measure.We consider now the duality result in the framework of the variable p in the space of Borel measures with L 2 (Ω) divergences.Surprisingly, the dual problem remains the same.We recall in this framework that as in (4.5), and d ≥ 1.Since we already assumed that α ∈ C(Ω) is positive, existence of a unique solution follows from Theorem 4.2.
We again propose to follow the Fenchel dual approach and let J : In this setting, it also holds that for every u ∈ L 2 (Ω), we have J * (div * u) = I K (u), i.e., Theorem 5.1.
In fact, we show that Lemma 5.3 and Lemma 5.4 remain true under the functional analytic setting of this section.
Proof of Lemma 5.3.
From Theorem 5.1, it follows that the duality result of Proposition 5.6 also holds in this setting; the proof is straightforward.

5.2.
The case when α is a measure.In this section, we will extend the duality result of Proposition 5.6 by letting α to be a non negative Borel measure, that is, α ∈ M + (Ω).However, p is a function in this setting.In its more general form, in problem (P * ), we set The results in this subsection are a generalization of the case of the Lebesgue integrable constraint α, that was presented in Section 5.1.We shall assume that p ∈ V Γ N (Ω), see (4.2) for the definition of V Γ N (Ω).
Since α ∈ M + (Ω), as we discussed in Section 1, we let the distributional gradient, and hence We prove that the dual problem to (P * ) is given by (P) with inequality constraint |Du| q ≤ α being understood in the sense of (2.5).
Recall that in section 3 we have shown existence and uniqueness of solution to (P).Here, we show that dual of problem (P * ) is given by (P).We start by writing (P * ) as (5.14) and F : L 2 (Ω) → R, as before, given by We prove now that Theorem 5.1 holds also true in the current setting.For brevity, we only discuss the essential modifications needed in Lemmas 5.3 and 5.4.
Proof of Theorem 5.1.This proof follows along the same lines as the proof to Theorem 5.1.We start by observing that the discussion in Remark 5.2 holds in the current setting as well, i.e., J * (div * u) only takes the values 0 and +∞.We now prove the result.The proof that J * (div * u) = 0 implies that u ∈ K follows along the lines of Lemma 5.3.Indeed (i) and (ii) in Lemma 5.3 apply directly, and for (iv) everything follows in the same way, when Du and dα are measures instead of the functions ∇u and α(x).
Finally, note that it follows identically as before that the polar function of F is given by Hence, the duality result of proposition 5.6 also holds in the case where α is a measure, with ∇ replaced by D.

A Finite Element Method with Applications
The purpose of this section is to illustrate the applicability of the proposed primal-dual approach to solve Problems (P) and (P * ).We assume throughout this section that p = q = 2.
Recall that Problem (P) in the case that α ∈ L ∞ (Ω) + is given by min and that the pre-dual problem (P * ) is given by min Now the first order (necessary and sufficient) optimality condition corresponding to (6.2) in the strong form is given by: Find where ∂ denotes the subdifferential operator.In order to solve (6.3), recall from the extremality conditions (5.13), that if u * and p * are solutions to (6.1) and (6.2), respectively, they satisfy u * := −div p + f a.e. in Ω .(6.4) Then, a primal-dual system arises from (6.3) and (6.4), which in the weak form becomes the following variational inequality of second kind: Due to their nonlinear and nonsmooth nature, it is challenging to solve (6.5)-(6.6).We shall proceed by introducing the Huber-regularization for φ(p) := |p| 2 in the last term under the integral in (6.2).This regularization is C 1 with piecewise differentiable first order derivative.Therefore one can use Newton type methods to solve the resulting regularized system.For a given parameter τ > 0, the Huber regularization of φ is given by As τ → 0, φ τ (p) → φ(p).Moreover, φ τ (•) is continuously differentiable with derivative given by 2) by φ τ (•), the regularized primal-dual system corresponding to (6.5)-(6.6) is given by (u, w) + (div p, w) = (f, w), for all w ∈ L 2 (Ω).( Notice, that φ τ (•) is piecewise differentiable and the second order derivative is given by where I d×d is the d × d identity matrix.
6.1.Finite Element Discretization.We discretize p and u using the lowest order Raviart-Thomas (RT 0 ) and piecewise constant (P 0 ) finite elements, respectively.Whenever needed, the integrals are computed using Gauss-quadrature which is exact for polynomials of degree less than equal to 4. For each fixed τ , we solve the discrete saddle point system (6.7)-(6.8)using Newton's method with backtracking line-search strategy.We stop the Newton iteration when each residual in L 2 (Ω)-norm is smaller than 10 −8 .Each linear solve during the Newton iteration is done using direct solve.Starting from τ = 10, a continuation strategy is applied where in each step we reduce τ by a factor of 1.30 until τ is less than equal to 10 −6 .We initialize the Newton's method by zero.To compute solution for next τ , we use the solution corresponding to previous τ as the initial iterate for the Newton's method.
6.2.Numerical Examples.Next, we report results from various numerical experiments.In all examples we consider Ω = (0, 1) × (0, 1) and we assume that Γ N = ∅, Γ D = ∂Ω, and hence pure Dirichlet boundary conditions on u on the entire boundary are set.In the first example, we construct exact solutions (p, u) when f and α are constants.We compare these exact solutions with our finite element approximation.These experiments validate our finite element implementation for constant α and f and provide optimal rate of convergence.Additionally, we solve (P * ) and (P) first for a fixed α and vary f and next we fix f and vary α.In our second experiment, we consider a more generic f with different features such as cone, valley and flat regions.In our final experiment, we consider α to be a measure.Notice that in this example, u is again Lipschitz continuous.In Figure 2 (top panel), we have shown the p − p h L 2 (Ω) and u − u h L 2 (Ω) when Ω = (0, 1) 2 , f = 1, and α = 1.We observe optimal rate of convergence in both cases.In the bottom row, the left panel shows u h , the middle panel shows |∇u h | 2 , and the right panel shows p h .We observe that, in this example, the gradient constraints are active in the entire region.Notice, that at the corners (which are sets of measure zero), gradient is undefined.  .Notice that we are touching the constraints in the entire region, except where the gradient is undefined.We have omitted the plots of the exact (u, p) as they look exactly same as (u h , p h ).
Next, we fix the number of p h and u h unknowns to be 197,120 and 131,072, respectively.Figure 3 shows our results for 3 different experiments.In all cases, we have used a fixed α = 1.The rows correspond to u h , |∇u h | 2 , and p h .Each column correspond to f = 1, f = 0.25, and f = 0.1.As expected, for a large value of f , we observe steep slope, but for smaller values of f , plateau regions appear.We also observe that the active region shrinks as f decreases since the gradient is zero at the top of the plateau.The dual variable p h also changes significantly with f .respectively.In all cases, we observe that the gradient constraints are active in the entire domain (except on a set of measure zero).For the case of piecewise constant α, nonsmoothness in p h is clearly visible.Notice that, the gradient constraints are active.Moreover, we also observe significant flat regions, where the gradient is zero.In Figure 7 have also displayed u h , |∇u h | 2 and p h when α = 1.5.The main novelty and challenge in this example is the fact that we let α to be a measure.Specifically ˆΩ v dα = ˆΩ v dx + 10 2 ˆω v dH 1 , for all v ∈ C ∞ c (Ω) and where ω := {(x, y) ∈ Ω : y = 0.5}, i.e., α consists in the Lebesgue measure dx and a weighted line measure on ω.
Let h denotes the meshsize, then α h is approximated as where ω h := {(x, y) ∈ Ω : 0.5 − 10 2 h ≤ y ≤ 0.5, x ∈ (0, 1)}.As h ↓ 0, we approximate the measure in the sense that ´Ω v dα h → ´Ω v dα for all v ∈ C ∞ c (Ω).When h = 8.4984 × 10 −5 , the results are shown in Figure 8 (top row).Finally, when h = 2.1412 × 10 −5 the results are provided in Figure 8 (bottom row).We notice that as h ↓ 0, we indeed approximate the measure: In fact, we observe a clear discontinuity on the solution u, the size of the jump is below 100 which is the upper bound on the distributional gradient on ω.We notice that as h ↓ 0, we accurately approximate the action of the measure α as a distributional gradient constraint.
is sequentially precompact in the sense of the intermediate convergence of BV(Ω) for any M > 0.Proof.Let {v n } be a sequence in K * , then it is bounded in BV(Ω), and thus v n v * in BV(Ω) for some v * ∈ BV(Ω) along a subsequence (not relabelled).Since |Dv n | p |Dv * | p in M(Ω), and |Dv n | p ≤ α it follows that for every open set O ⊂ Ω that |Dv * | p (O) ≤ lim inf n→∞ |Dv n | p (O) ≤ α(O), (3.1)where we have used the lower-semicontinuity property for open sets of weak convergence of measures; see [12, Proposition 4.2.3].Additionally, since elements in M(Ω) are outer (and inner) regular, we have that for a Borel set B it holds that µ(B) = inf µ(O) where the infimum is taken over all open sets such that B ⊂ O; see [12, Proposition 4.2.1].Thus, |Dv * | p (B) ≤ α(B) (3.2) follows from (3.1) by taking the infimum over {O open : B ⊂ O}.In order to prove that v n converges to v * in the sense of intermediate convergence, we are only left to prove that |Dv n | |Dv * | narrowly in M + (Ω) (see [12, Proposition 10.1.2]).The latter meaning that ´Ω ϕ|Dv n | → ´Ω ϕ|Dv * | for each continuous and bounded ϕ on Ω.Given that α ∈ M + (Ω) we have that for each > 0 there exists a compact set Λ ⊂ Ω such thatα(Ω \ Λ ) ≤ .Since v n ∈ K, then |Dv n | ≤ α,and hence for each > 0 the compact set Λ ⊂ Ω, is such that |Dv n |(Ω \ Λ ) ≤ , for all n ∈ N.Then, by Prokhorov Theorem (see [12, Theorem 4.2.3]),there is a subsequence of {|Dv n |} (not relabelled) that |Dv n | |Dv * | narrowly in M + (Ω).That is, along a subsequence, v n converges to v * in the sense of intermediate convergence.This implies that v * ∈ BV Γ D (Ω), by virtue of Theorem 2.3 and the fact that v n ∈ BV Γ D (Ω) for all n ∈ N.

Figure 2 . 1 :
Figure 2. Example 1: Top panel -We have shown the L 2 (Ω)-error between the computed solution (u h , p h ) and the exact solution (u, p).The optimal linear rate of convergence is observed.Bottom panel: Computed u h (left), |∇u h | 2 (middle), p h (right).Notice that we are touching the constraints in the entire region, except where the gradient is undefined.We have omitted the plots of the exact (u, p) as they look exactly same as (u h , p h ).

Figure 3 .
Figure 3. Example 1 (fixed α, varying f ).The rows correspond to u h , |∇u h | 2 , and p h .The columns represent f = 1, 0.25, and 0.1.In all cases, we observe that the gradient constraints are active but the activity region shrinks as f decreases, this is expected since the gradient on the plateau region is zero.The behavior of p also changes considerably with f .In Figure 4 we again show results from 3 different experiments.In all cases, we have used a fixed f = 1.The rows correspond to u h , |∇u h | 2 , and p h .Each column corresponds to α = 1, α = 0.5, and α = 0.75, y ≤ 1 − x, 1.0, otherwise,

Figure 4 .Example 2 . 2 +
Figure 4. Example 1 (fixed f , varying α).The rows correspond to u h , |∇u h | 2 , and p h .The first two columns represent constant α = 1 and 0.5.The third column corresponds to α with jump discontinuity.In all cases, we observe that the gradient constraints are active in the entire region, except on a set of measure zero.Moreover, discontinuity in p in the last column is clearly visible.

Figure 8 .
Figure 8. Example 4 (α measure).Top row: u h and p h when the meshsize h = 8.4984 × 10 −5 .Bottom row: u h and p h when the meshsize h = 2.1412 × 10 −5 .We notice that as h ↓ 0, we accurately approximate the action of the measure α as a distributional gradient constraint.
Proposition 2.4.The condition in (2.2) is equivalent to ˆΩ w|Dv| p ≤ ˆΩ w dα for every w ∈ C ∞ (Ω) with w ≥ 0 in Ω. (2.5) Proof.Suppose that (2.2) holds true and let K n be a sequence of closed sets such that ˆΩ\Kn |Dv| p → 0 and The sequence {K n } exists given that measures in M + (Ω) are inner regular; see [12, Proposition 4.2.1].Let w ∈ C ∞ (Ω) be nonnegative and arbitrary.Accordingly, let {w n } in C ∞ 0 (Ω) be nonnegative, uniformly bounded in Ω, and such that w n = w in K n .Hence | w| + |w n | can be uniformly estimated by a constant, and by (2.6) it holds that ˆΩ( w − w n 1)in(1.3)is understood in the sense of measures, i.e., (2.1) holds true if ˆΩ w|Dv| p ≤ ˆΩ w dα for all w ∈ C ∞ for all open sets O ⊂ Ω.It is possible to replace C ∞ 0 (Ω) in (2.2) by C ∞ (Ω), which we discuss next.