The pursuit of a dream, Francisco Javier Sayas and the HDG methods

Franciso Javier Sayas, man of grit and determination, left his hometown of Zaragoza in 2007 in pursuit of a dream, to become a scholar in the USA. I hosted him in Minneapolis, where he spent three long years of an arduous transition before obtaining a permanent position at the University of Delaware. There, he enthusiastically worked on the unfolding of his dream until his life was tragically cut short by cancer, at only 50. In this paper, I try to bring to light the part of his academic life he shared with me. As we both worked on hybridizable discontinuous Galerkin methods, and he wrote a book on the subject, I will tell Javier’s life as it developed around this topic. First, I will show how the ideas of static condensation and hybridization, proposed back in the mid 60s, lead to the introduction of those methods. This background material will allow me to tell the story of the evolution of the hybridizable discontinuous Galerkin methods and describe Javier’s participation in it. Javier faced death with open eyes and poised dignity. I will end with a poem he liked.


Javier's Minneapolis triptych
The second triptych is about Javier's life in Minneapolis.
The top panel depicts El descubrimiento de América por Cristóbal Colón, oil on canvas by Salvador Dalí, 1958/59. The many crosses are Javier's Western, Christian cultural forces which propelled him across the Atlantic. And Gala stands for the promise of tenure in the American New World. Javier arrives to the New Continent with the anticipation of the encounter with an academic El Dorado and a trepidation, for stepping into uncharted territory, somewhat assuaged by a little help from his friends: Gabriel Gatica, from the University of Concepción, and Tomás Chacón, from the University of Sevilla, had alerted me about Javier's transatlantic venture. Harry Singh, from the School of Mathematics of the University of Minnesota, Rosario Grau (my wife) and I were ready to host him in Minneapolis (Fig. 2).
The panel of the middle is Tránsito en espiral, oil on canvas by Remedios Varo, 1962. Remarkable woman of strong political convictions, Remedios Varo was a student of G.I. Gurdieff's ideas and one of the main surrealist painters of her time. Like Javier, she was Spaniard, immigrated to the New World (México), enjoyed 9 years of success as an artist and also died tragically, at 54.
Remedios Varo could very well understand all the works and tribulations that assail an immigrant. Few have seen the actual Wall, the energy barrier to the completion of a transitional but long, sustained and difficult task, as well as Remedios Varo did. As we see in this panel, from his basis in Minneapolis, Javier put three long years to circumnavigate the Wall with great effort, grit and determination. Note how Javier is pedaling furiously strange-looking vehicles along the channel made by the spiraling wall. Some of them are papers by himself, like the 2009 SINUM paper [59] whose main result was celebrated by re-publishing the paper later in 2016 in the prestigious SIAM Review [60]. Others are papers with his new colaborators Jay Gopalakrishnan, Johnny Guzmán, Youngmock Jeon, Ngoc-Cuong Nguyen, Jaime Peraire and Manuel Solano. Yet others, those smaller and closer to the tower of the center, are his no less relevant job interview talks.
The panel of the bottom shows Saint George defeating the Dragon, by Johann König, 1630. It signals that the spiraling iterations around the tower are over and that the desired goal is finally achieved. We see Javier giving a great interview talk at the University of Delaware thus defeating the dark forces that had been barring his entrance to the world of scholars in the USA. On the left and behind him, a crowned tenure-track is waiting for him. At this point, Javier knew that his dream would be realized (Fig. 3).

Javier's Delaware triptych
The third and last triptych has the classical configuration and displays the covers of the three books that Javier wrote and published during his tenure at the University of Delaware. They are a reflection of Javier's passion for the craft of clarification and mathematical exposition, and represent the rich academic life that Javier was able to realize.
The book on the right, written with his then Ph.D. student Shukai Du, is about the socalled hybridizable discontinuous Galerkin (HDG) methods [38]. It was completed only a few weeks before Javier's death. To honor Javier's interest in the HDG methods, a favorite topic of mine, in what follows, I describe these methods and tell the story of how Javier participated in their evolution.

Static condensation and hybridization: devising HDG methods
We present the HDG methods, not as they were originally introduced in 2009 [20], but by closely following the more didactic narrative which can be found in the 2016 review [6], as it merges the development of these methods with the evolution of the ideas of static condensation [47] and hybridization [42] introduced in 1965.
These two ideas were proposed in the framework of numerical methods for linear elasticity problems and were originally thought of as techniques for improving the matrices associated to the weak formulations of the methods, that is, for making them smaller and easier to numerically invert. Next, we show that we can also view those ideas as a single property of the exact solution which can then be used as a template to be applied to any of its approximations, like the continuous Galerkin and the mixed methods. Such template will be used for devising the HDG methods.

The model problem and its exact solution
We proceed in the framework of steady-state diffusion problems. So, we consider the following second-order elliptic model problem: Here c is a matrix-valued function which is symmetric and uniformly positive definite on , f is an L 2 ( ) function, and · represents the trace operator. For any element K of a mesh T h of the domain , we know that the exact solution satisfies the local problem We also know that for any interior face F = ∂ K + ∩ ∂ K − , K + , K − ∈ T h , it also satisfies the transmission conditions Here, η ± represents the trace, of the generic function η, on the face F from the interior of the element K ± , and n ± represents the unit normal to F pointing outward from K ± . Note that the first transmission condition is a reflection of the presence of the term ∇u in the first equation. The second transmission condition is a reflection of the presence of the term ∇ · q in the second equation. Finally, we trivially know that the Dirichlet boundary condition is also satisfied.
Conversely, suppose that we are given a sufficiently smooth function u defined on F h , the set of all faces F of the mesh T h . If for each element K of the mesh T h , we define (q, u) in terms of u and f as the solution of the local problem we can then ask ourselves for what funcion u we can say that (q, u) is the actual exact solution. The answer is that the function u must be the solution of the global problem where q is the trace of q = q( u, f ) on ∂ K . This is the property of the exact solution we alluded to above. Expressing (q, u) in terms of the unknown u (and f ) corresponds to the static condensation and the introduction of the unknown u corresponds to the hybridization. We can further refine this property by separating the influence of u from that of f . Indeed, by the principle of superposition, we have that (q, and u is the solution of the equations Example Let us illustrate this characterization of the exact solution in the one-dimensional and where the function u is the solution of the global problem A simple computation gives us the solution of the local problem on (x i−1 , x i ), where G i (·, ·) is the Green's function of the local problem (see §2.1.2 in [6]), G i x (·, ·) its partial derivative with respect to the first variable, and ϕ i is the piecewise-linear function that is equal to one at x i and zero at x j for j = i. The function u is then the solution of Note that u i is nothing but the value of the exact solution at the partition point x i . So the characterization of the exact solution shows that, if we can solve all the local problems exactly, the remaining values of the exact solution, u, can be obtained by solving a matrix equation.

The continuous Galerkin method
Now, let us apply the above template to the well-known continuous Galerkin method. This method provides an approximation to u, solution of our model problem rewritten as where a := c −1 . Now, let us set where, for each element K ∈ T h , W (K ) is a finite dimensional space (typically of polynomials), and I h is an interpolator or projection operator into the space of approximate traces {w| : w ∈ W h }. Then, we take the approximation u h as the element of W h (u D ) determined as the solution of the weak formulation To apply the template, for each element K ∈ T h , we split the local space W (K ) as follows: This induces the following split in the global spaces: and allows us to rewrite the weak formulation defining u h as follows. First, we obtain U ∈ W (K ) in terms of u h and f by solving Then we have that U = u h if and only if we take the function and define it to be the solution of The fact that the first of the equations determining u h can be interpreted as a transmission condition is proved in detail in [6]. There, we can also find the following further refinement of the above result. We have that [47]. What we just did by using weak formulations, Guyan [47] did with matrices back in 1965. Let us establish a parallel between these two versions of the same procedure. The system of equations associated to the original weak formulation of the continuous Galerkin method is

Relation with the static condensation technique by Guyan
After splitting the degrees of freedom, we can write that

Then the original equation becomes
The matrix form of the solution of the local problems is . and the matrix form of the transmission condition is Note that the process of static condensation just described in terms of matrices is nothing but a block-matrix Gauss elimination, and that the matrix associated to the transmission condition is nothing but the corresponding Schur complement matrix.
Note also that the unknown u h is completely captured by the part of the approximate solution we called U ∂ . This function already existed as part of the approximation and so the process of hybridization remained unnoticed. Indeed, no new, hybrid unknown had to be introduced. This is why only static condensation is attributed to Guyan [47]. The process of hybridization became evident in the framework of mixed methods, as we show next.
Example But before doing that, let us consider the application of the template to the one-dimensional case. For W (K ) := P k (K ), the solutions of the local problems are is the discrete Green's function (see §2.2.4 in [6]), and the global problem for the values are exact in this case as we can see if we compare this equation with the corresponding one for the exact solution.

A mixed method defines an approximation to
The scalar-valued functions in W h can be discontinuous across inter-element boundaries. However, the vector-valued functions in V h must have continuous normal components across inter-element boundaries. We assume that, for every element K ∈ T h , the finite dimensional local spaces V (K ) and W (K ) (typically of polynomials) have been chosen so that the method is well defined.
To apply the template under consideration to this method, we first have to hybridize the method by introducing a new unkown u h which we take in a suitably chosen space M h . Only then we will be able to statically condense the method. To show how to do that, we follow [17].
We define (Q, U) ∈ V (K ) × W (K ) in terms of u h and f as the solution of the local problem Then, we take u h as the element of the space M h (u D ) satisfying Note that we can impose the Dirichlet boundary condition strongly since we do not take the restriction of M h to ∂ to lie in a finite dimensional space.
It is not difficult to show that this choice of u h is the only one which ensures that (Q, U) = (q h , u h ), see [17] for details. Therein, it has also been proved that we can further refine the above characterization of the mixed method approximate solution as follows. We have that and the function u h is the element of M h (u D ) which solves the global problem

Relation with the hybridization of Fraejis de Veubeke [42]
Let where [ u h ] is the vector of the degrees of freedom of the hybrid unknown on the interior faces. This system of equations is significantly bigger than the original one, but is such that Note, in contrast, the symplicity of the corresponding weak formulation! Once again, we can see that the process of static condensation just described in terms of matrices is nothing but a block-matrix Gauss elimination, and that the matrix associated to the transmission condition is nothing but the corresponding Schur complement matrix.
Note also that the unknown u h was not part of the approximate solution of the mixed method. Its introduction marks the effective application of what is called hybridization. At first, it increases significantly the size of the matrix equation, but then it allows for an efficient static condensation as it results in a smaller (and symmetric, positive definite!) matrix for the hybrid unknown.
We are now ready to apply the template for the devising of the HDG methods. But before, let us quickly see how this process looks in the simple one dimensional case we have been considering.
Example Let us now show the mixed method in the one-dimensional case. For the wellknown Raviart-Thomas (RT) [56] space of index k ≥ 1, namely, for V (K ) × W (K ) := P k+1 (K ) × P k (K ), the solutions of the local problems are where G i h (·, ·) is the discrete local Green's function given by the method and −H i h (·, ·) its corresponding discrete x-derivative (see §2. 3

.3 in [6]), and the global problem for the values
Again the values { u i } N i=0 are exact.

HDG methods
We now use the template to devise the HDG methods. The main advantage of this approach is that the static condensation and hybridization of the resulting methods are guaranteed by construction, see [20].
As we see next, the distinctive features of the HDG methods are that (i) the local problems are expressed in terms of discontinuous Galerkin methods, and that (ii) the transmission condition characterizing the hybrid unkown u h are nothing but the requirement (standard for discontinuous Galerkin methods) that the numerical flux be conservative.
So, on the element K ∈ T h , we define the approximation given by an HDG method, (q h , u h ), in terms of ( u h , f ) as the element of V (K )× W (K ) which solves the discontinuous Galerkin (DG) scheme (3.1) The hybrid unknown u h is defined as follows. For each face F ∈ F i h , we take u h | F in a suitably defined finite dimensional space M(F). Then, we determine u h as the solution of This completes the devising of the HDG methods by using this first template of the exact solution.
Note that, as wanted, the only globally-coupled degrees of freedom are those of the hybrid unknown, u h . Note also that, for each interior face F, there are two values of the numerical flux q h . The transmission condition, that is, the first equation determining the hybrid unknown u h , forces the normal components of those two values to weakly have the same value. In this way, a weak version of the local conservativity property, typical of DG methods, is imposed.
Suppose now that, on each face of an element, the stabilization function τ is the simple multiplication by the parameter τ . Then, if the transmission condition implies that [[ q h ]] = 0, we have that provided τ + + τ − > 0. These numerical traces are a particular case of the numerical traces of the DG methods considered in [3,35]. Those DG methods are hybridizable and hence amenable to an efficient static condensation. This is what gave the name HDG to these methods.
Note also that the hybrid function u h is both data for the local problems as well as the unknown for the global problem. The hybridization and static condensation approach used here to devise HDG methods for the steady-state diffusion problem can be extended to any partial differential equation by using what we could call a generalized hybridization and static condensation procedure. Indeed, the hybrid unknowns are nothing but any data which would render the local problems well posed. The global problem would then be the one that determines the hybrid unknowns. See the 2005 review [18] on hybridization techniques, and the 2009 paper [19] where HDG methods for the Stokes equations are devised by using four different choices of the hybrid unknowns.
As we did for the continuous Galerkin and the mixed methods, we can further refine the presentation of the HDG methods and rewrite them by separating the influence of u h from that of f . We can express the approximation (q h , u h ) as the sum and the numerical trace u h is the element of the affine space satisfying the equations Here, (·, ·) T h := K ∈T h (·, ·) K and ·, · ∂T h := K ∈T h ·, · ∂ K . We can then see that the matrix equation for u h is symmetric. That it is also positive definite is easy to show when we assume that For a proof, see [6,20]. In sharp contrast with the mixed methods, note how simple is to choose the local spaces of the HDG methods. The HDG methods inherit from the DG methods their ability of handling elements of arbitrary shape and arbitrary basis functions, as no strong inter-element continuity is required. The different HDG methods associated to this template, are generated by choosing the local spaces V (K ), W (K ), M(F) and the stabilization function τ . Accordingly, two lines of research were immeditely generated which are still producing interesting results. One of them focuses on how to choose the local spaces whereas the other focuses on how to chose the stabilization functions.

Superconvergent HDG methods
Most of the work that Javier did on the HDG methods was on the line of research which emphasizes the choice of the local spaces. This approach takes advantage of the similarity of the formulations of the hybridized version of the mixed methods and the HDG methods. Indeed, note that, if in the latter formulation we set the stabilization function τ to zero, we do get a hybridized mixed method. This is why a natural question to ask is what are the local spaces for which an HDG method has the same convergence properties of a mixed method? In particular, it is very well known that the RT method [56] for simplexes is such that a projection of the error u − u h converges with an additional power of h. Using this, one can devise an element-by-element post processing leading to a new approximation u * h which converges with an additional power of h. That is why one says that this mixed method superconverges. The question now is if there are HDG methods which superconverge?
The first affirmative answer was given in 2008 [8] where the single face hybridizable (SFH) method was introduced. Defined on simplexes, the stabilization function was set to be zero on all faces of the simplex except on one, arbitrarily chosen. The second superconvergent HDG methods were found also in 2009 [24]. They use local discontinuous Galerkin (LDG) methods to define the local problems, are applied on meshes made of simplexes and use stabilization functions uniformily bounded above and below. This result was obtained by exploiting the above-mentioned similarity between the formulations of the HDG and the mixed methods.
In the 2010 paper [22], of which Javier was one of the authors, the first projection-based analysis of superconvergent HDG methods, defined on simplexes, was obtained. Such type of analysis was prompted by similar ones for the mixed methods. The mixed method projections ensure a powerful commuting diagram property which was thought to be necessary for their superconvergence property. The result in [22] showed that superconvergence can take place without such commuting diagram property. From this point on, all analyses of HDG methods would be projection-based.
Javier also participated in the extention of this projection-based technique for the analysis of superconvergent HDG methods for the Stokes flow of incompressible flow in 2011 [21] and again in 2014 [30] where the much more difficult case of exactly divergence-free approximate velocities was treated. The projection-based analysis was also extended for HDG methods for Timoshenko beams in 2012 [4].
All the above superconvergent HDG methods were obtained for simplicial elements. The question was now if there were superconvergent HDG methods defined in elements of arbitrary shape. The search for sufficient conditions for an HDG method to be superconvergent was carried out in 2012 [26,27], for steady-state diffusion problems. The extension of the projection-based analysis proposed in [22] was caried out to all these methods, and new superconvergent HDG methods for squares, cubes and prisms were constructed. The extension of this result to steady-state linear elasticity was done in 2013 in [33], and to the Stokes equations of incompressible flow also in 2013 in [32]; see the 2014 review [34].
At this point, no superconvergent HDG method for quadrilaterals or general, flat-faced hexahedra had been constructed. The work we just mentioned on the sufficient conditions for superconvergence had to be slighly tightened up for that to take place. So, in 2017, a theory of M-decompositions was introduced [14]. Javier did participate in its creation. Its application to the construction of superconvergent HDG methods defined on general polygonal elements was carried out in 2017 in [9] and for superconvergent HDG methods defined on general pyramids, prisms, and hexahedra also in 2017 in [10]. The extension to Stokes flow was done in 2017 in [13], and to elasticity in 2018 in [12]. The theory of M-decompositions can deliver new superconvergent mixed methods as it can construct superconvergent HDG methods. So, it can also be considered as a technique for constructing mixed methods for elements of arbitrary shape. As such, it can be used for the systematic construction of commuting exact sequences on elements of the above-mentioned shapes, see the 2017 paper [11]. For an introductory overview of the theory of M-decompositions, see [15].

Coupling HDG and boundary element methods
During his Zaragoza period, Javier had been interested in the coupling of LDG and BE methods [2,45,46]. During his Minneapolis period, he continued to develop this interest in three papers published in 2012, namely, on the coupling of the RT or HDG with BE methods [23], on general symmetric couplings of DG and BE methods [29], and on coupling at a distance HDG and BE methods [31]. During his Delaware period, he worked on nonsymmetric couplings of HDG and BE methods in the 2017 paper [44].

Wave propagation
The first HDG methods for the acoustic and elastic wave equations were implicit and were proposed in 2011 [51]. In 2014, the semidiscrete version of the methods for the acoustic wave equation was proven to be uniformly-in-time superconvergent in [28]. In 2016, explicit HDG methods for wave propagation were considered in [61]. Also on that year, an overview of the state of art on HDG method for hyperbolic problems appeared [25] which describes the state of the art up to that time.
It is, more or less, around this time that Javier took an interest in HDG methods for wave propagation. Javier developed MATLAB software for the HDG method in 3D, see the 2015 paper [43], which was used, for example, in his 2018 paper on viscoelastic waves [1]. During that period, he also became interested in the application of HDG methods to wave propagation phenomena. Javier worked on the Stormer-Numerov HDG methods for the wave equation which was published in 2018 [16]. This was the last paper on which we collaborated. The results of the paper prompted the introduction of the Symplectic-Hamiltonian HDG methods: First for the acoustic wave equation, in 2017 [57], and then for the equations of elastodynamics, in 2021 [58]. Javier developed HDG methods for the time-harmonic linear elasticity in 2017 [48].

Unified analyses, new stabilizations
Javier also did work on the other main line of research on HDG methods, namely, the search for stabilization functions. In [49], see also [50], what is called the Lehrenfeld-Schöberl projection was introduced. It allows to define new HDG methods converging optimally regardless of the shape of the elements. The LDG-H methods also enjoy this property, but the spaces M h allowed with the use of the above-mentioned projection are smaller. In 2015 [54] we can find an analysis of the corresponding methods for steady-state diffusion; see also [55] the case of Stokes flow.
Another contribution in this line of research is the introduction of the hybrid high-oder (HHO) methods in 2014 [36,37]. For these methods, optimal convergence and superconvergence properties are achieved for elements of general shape. In 2016 [7] it was shown that the HHO method can be rewritten as an HDG method with a special stabilization. A distinctive role of this stabilization is to provide a higher-order of accuracy which the standard stabilizations do not provide.
Javier became interested in unifying the analyses of the HDG methods using the old stabilizations and the HDG using the Lehrenfled-Schöberl stabilization, which he called the HDG+ methods. He worked on a unified analysis of HDG methods for the steady-state Maxwell's equations in a paper published in 2020 [40]; notably, it included the HDG+ method introduced in 2017 in [5]. Also published in 2020, the paper [39] presents a single, simple and concise analysis of several variants of the HDG+ method using a projection-based approach. An interesting novelty here is that the theory of M-decompositions is used, but only as a intermediate step. Javier thus brought together the two lines of research on HDG methods. Finally, published in 2021 is a paper [41] in which a simple way of constructing HDG+ projections on polyhedral elements is introduced. This was Javier's last paper on HDG methods.

Epilogue: the little poem of things that statically condense
Javier's legacy not only lies in the many papers he wrote and the many Ph.D. Students he graduated. He touched many of us by his willingness to be actively engaged with the mathematical community. The yearly DelMar (Delaware-Maryland) Numerics Seminar started in 2012 and was originally organized by Ricardo H. Nochetto (University of Maryland), Petr Plechac and Francisco-Javier Sayas (University of Delaware), and by Tobias von Petersdorff (University of Maryland). It is thus fitting that, in memory of Javier, this seminar is now called the Sayas Numerical Day.
Moreover, the Sayas Numerics Seminar is a weekly online seminar which emphasizes, as Javier would have wanted, the participation of early career researchers. It has run during Fig. 4 The little poem of things that statically condense. The poem is traditionaly, but wrongly, attributed to Moments of the Soul, Monastery of Bgheno-Noravank in Syunik, 989. Javier's picture is taken from [52] What does statically condense? When implementing HDG methods, the degrees of freedom, When drinking in a seminar BBQ, dew-like drops on a beer, When saying farewell to a friend, tears on the heart.
the Fall semester of 2020 and the Spring semester of 2021, and I foresee it running many more semesters to come. Its organizers are Harbir Antil (George Mason University), Andrei Draganescu (University of Maryland), Ricardo H. Nochetto, Petr Plechac and Tobias von Petersdorff. The picture in Fig. 4 is the one that appears on the corresponding webpage. Its creation was instigated by Ricardo Nochetto [53] and carried out by Tobias von Petersdorff [62] by taking Javier's picture from the obituary on the University of Delaware website and then applying on it an adaptivity Matlab code he wrote. To me, Javier was first a guest, then a collaborator and finally a friend. I want to leave in writing, verba volant, scripta manent, what I said to him on Saturday February 16th, 2019, in his meeting Numerical Analysis & No Regrets: Farewell Javier! It was good to play with you.
To end, let me summarize in Fig. 4 the story told in this short paper with a poem for Javier I read on that Saturday: The little poem of things that statically condense. We joked about the attribution. He then told me he liked the poem, accused me of being emotional, and gave me a big hug.
Chairman, for the invitation to participate in the XXVI CEDYA/ XVI CMA, June 14-18, 2021. CEDYA stands for Congreso de Ecuaciones Diferenciales y Aplicaciones, and CMA for Congreso de Matemática Aplicada. Concerning the paper itself, let me thank the referees for their constructive comments, and Gabriel Gatica for completing the information I had on the top panel of the first of Javier's triptychs. I would also like to thank Ricardo Nochetto and Tobias von Petersdorff for letting me use their figure of Javier [52] which I am displaying in Fig. 4. Thanks are due too to Scot Adams, Yanlai Chen, Shukai Du and Salim Meddahi for their sharp proofreading, and then again to Shukai Du for help in properly citing some of his work with Javier. Last but not least, I would like to thank Rosario Grau for her lucid editorial remarks.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.