Leszek Demkowicz

Ali Vaziri Astaneh, Federico Fuentes, Jaime Mora et al.

This work represents the first endeavor in using ultraweak formulations to implement high-order polygonal finite element methods via the discontinuous Petrov-Galerkin (DPG) methodology. Ultraweak variational formulations are nonstandard in that all the weight of the derivatives lies in the test space, while most of the trial space can be chosen as copies of $L^2$-discretizations that have no need to be continuous across adjacent elements. Additionally, the test spaces are broken along the mesh interfaces. This allows one to construct conforming polygonal finite element methods, termed here as PolyDPG methods, by defining most spaces by restriction of a bounding triangle or box to the polygonal element. The only variables that require nontrivial compatibility across elements are the so-called interface or skeleton variables, which can be defined directly on the element boundaries. Unlike other high-order polygonal methods, PolyDPG methods do not require ad hoc stabilization terms thanks to the crafted stability of the DPG methodology. A proof of convergence of the form $h^p$ is provided and corroborated through several illustrative numerical examples. These include polygonal meshes with $n$-sided convex elements and with highly distorted concave elements, as well as the modeling of discontinuous material properties along an arbitrary interface that cuts a uniform grid. Since PolyDPG methods have a natural a posteriori error estimator a polygonal adaptive strategy is developed and compared to standard adaptivity schemes based on constrained hanging nodes. This work is also accompanied by an open-source $\texttt{PolyDPG}$ software supporting polygonal and conventional elements.

Weifeng Qiu, Leszek Demkowicz

The paper presents a generalization of Arnold-Falk-Winther elements for three dimensional linear elasticity, to meshes with elements of variable order. The generalization is straightforward but the stability analysis involves a non-trivial modification of involved interpolation operators. The analysis addresses only the h-convergence.

Jamie Bramwell, Leszek Demkowicz, Jay Gopalakrishnan et al.

We present two new methods for linear elasticity with simultaneously yield stress and displacement approximations of optimal accuracy in both the mesh size h and polynomial degree p. This is achieved within the recently developed discontinuous Petrov-Galerkin (DPG) framework. In this framework, both the stress and the displacement approximations are discontinuous across element interfaces. We study locking-free convergence properties and the interrelationships between the two DPG methods.

Brendan Keith, Federico Fuentes, Leszek Demkowicz

The flexibility of the DPG methodology is exposed by solving the linear elasticity equations under different variational formulations, including some with non-symmetric functional settings (different infinite-dimensional trial and test spaces). The family of formulations presented are proved to be mutually ill or well-posed when using traditional energy spaces on the whole domain. Moreover, they are shown to remain well-posed when using broken energy spaces and interface variables. Four variational formulations are solved in 3D using the DPG methodology. Numerical evidence is given for both smooth and singular solutions and the expected convergence rates are observed.

Brendan Keith, Socratis Petrides, Federico Fuentes et al.

A finite element methodology for large classes of variational boundary value problems is defined which involves discretizing two linear operators: (1) the differential operator defining the spatial boundary value problem; and (2) a Riesz map on the test space. The resulting linear system is overdetermined. Two different approaches for solving the system are suggested (although others are discussed): (1) solving the associated normal equation with linear solvers for symmetric positive-definite systems (e.g. Cholesky factorization); and (2) solving the overdetermined system with orthogonalization algorithms (e.g. QR factorization). The finite element assembly algorithm for each of these approaches is described in detail. The normal equation approach is usually faster for direct solvers and requires less storage. The second approach reduces the condition number of the system by a power of two and is less sensitive to round-off error. The rectangular stiffness matrix of second approach is demonstrated to have condition number $\mathcal{O}(h^{-1})$ for a variety of formulations of Poisson's equation. The stiffness matrix from the normal equation approach is demonstrated to be related to the monolithic stiffness matrices of least-squares finite element methods and it is proved that the two are identical in some cases. An example with Poisson's equation indicates that the solutions of these two different linear systems can be nearly indistinguishable (if round-off error is not an issue) and rapidly converge to each other. The orthogonalization approach is suggested to be beneficial for problems which induce poorly conditioned linear systems. Experiments with Poisson's equation in single-precision arithmetic as well as the linear acoustics problem near resonance in double-precision arithmetic verify this conclusion.

Ali Vaziri Astaneh, Brendan Keith, Leszek Demkowicz

In this article, several discontinuous Petrov-Galerkin (DPG) methods with perfectly matched layers (PMLs) are derived along with their quasi-optimal graph test norms. Ultimately, two different complex coordinate stretching strategies are considered in these derivations. Unlike with classical formulations used by Bubnov-Galerkin methods, with so-called ultraweak variational formulations, these two strategies in fact deliver different formulations in the PML region. One of the strategies, which is argued to be more physically natural, is employed for numerically solving two- and three-dimensional time-harmonic acoustic, elastic, and electromagnetic wave propagation problems, defined in unbounded domains. Through these numerical experiments, efficacy of the new DPG methods with PMLs is verified.

Jaime Mora, Leszek Demkowicz

Numerical integration of the stiffness matrix in higher order finite element (FE) methods is recognized as one of the heaviest computational tasks in a FE solver. The problem becomes even more relevant when computing the Gram matrix in the algorithm of the Discontinuous Petrov Galerkin (DPG) FE methodology. Making use of 3D tensor-product shape functions, and the concept of sum factorization, known from standard high order FE and spectral methods, here we take advantage of this idea for the entire exact sequence of FE spaces defined on the hexahedron. The key piece to the presented algorithms is the exact sequence for the one-dimensional element, and use of hierarchical shape functions. Consistent with existing results, the presented algorithms for the integration of $H^1$, $H(\text{curl})$, $H(\text{div})$, and $L^2$ inner products, have the $O(p^7)$ computational complexity. Additionally, a modified version of the algorithms is proposed when the element map can be simplified, resulting in the reduced $O(p^6)$ complexity. Use of Legendre polynomials for shape functions is critical in this implementation. Computational experiments performed with $H^1$, $H(\text{div})$ and $H(\text{curl})$ test shape functions show good correspondence with the expected rates.

Brendan Keith, Leszek Demkowicz, Jay Gopalakrishnan

We introduce a cousin of the DPG method - the DPG* method - discuss their relationship and compare the two methods through numerical experiments.

Weifeng Qiu, Leszek Demkowicz

We continue our study on variable order Arnold-Falk-Winther elements for 2D elasticity in context of both affine and parametric curvilinear elements. We present an $h$-stability result for affine elements, and an asymptotic stability result for curvilinear elements. Both theoretical results are confirmed with numerical experiments.

Jaime Mora-Paz, Stefan Henneking, Leszek Demkowicz et al.

With the goal of accurately extracting the optical field losses in a three-dimensional (3D), circularly coiled waveguide (e.g., bent optical fiber), this effort presents the numerical methodologies that are implemented for an envelope Maxwell model that propagates electromagnetic fields as an entirely boundary value problem. Our unique modeling approach includes an ultraweak variational formulation of the envelope Maxwell model in the curved geometry of the bending, which is discretized by the discontinuous Petrov-Galerkin (DPG) method, which permits residual-driven mesh and polynomial-order adaptivity. This also, then, requires a unique approach for constructing perfectly matched layers (PMLs) as absorbing boundary conditions in both the direction of optical field propagation and in the tangential directions, where unguided energy escapes the waveguide. Our coiled waveguide modeling technology extracts the mode confinement losses from the propagation of the coherent optical field through the bent waveguide. We verify our simulations against the semi-analytical results from the analogous bent slab waveguide problem, and we successfully demonstrate stable convergence to loss values for the 3D coiled optical fiber problem, which has never been done previously for our specific modeling approach.

Kamil Doległo, Anna Paszyńska, Maciej Paszyński et al.

This paper deals with the following important research question. Traditionally, the neural network employs non-linear activation functions concatenated with linear operators to approximate a given physical phenomenon. They "fill the space" with the concatenations of the activation functions and linear operators and adjust their coefficients to approximate the physical phenomena. We claim that it is better to "fill the space" with linear combinations of smooth higher-order B-splines base functions as employed by isogeometric analysis and utilize the neural networks to adjust the coefficients of linear combinations. In other words, the possibilities of using neural networks for approximating the B-spline base functions' coefficients and by approximating the solution directly are evaluated. Solving differential equations with neural networks has been proposed by Maziar Raissi et al. in 2017 by introducing Physics-informed Neural Networks (PINN), which naturally encode underlying physical laws as prior information. Approximation of coefficients using a function as an input leverages the well-known capability of neural networks being universal function approximators. In essence, in the PINN approach the network approximates the value of the given field at a given point. We present an alternative approach, where the physcial quantity is approximated as a linear combination of smooth B-spline basis functions, and the neural network approximates the coefficients of B-splines. This research compares results from the DNN approximating the coefficients of the linear combination of B-spline basis functions, with the DNN approximating the solution directly. We show that our approach is cheaper and more accurate when approximating smooth physical fields.

Leszek Demkowicz, Jay Gopalakrishnan, Brendan Keith

This article introduces the DPG-star (from now on, denoted DPG$^*$) finite element method. It is a method that is in some sense dual to the discontinuous Petrov-Galerkin (DPG) method. The DPG methodology can be viewed as a means to solve an overdetermined discretization of a boundary value problem. In the same vein, the DPG$^*$ methodology is a means to solve an underdetermined discretization. These two viewpoints are developed by embedding the same operator equation into two different saddle-point problems. The analyses of the two problems have many common elements. Comparison to other methods in the literature round out the newly garnered perspective. Notably, DPG$^*$ and DPG methods can be seen as generalizations of $\mathcal{L}\mathcal{L}^\ast$ and least-squares methods, respectively. A priori error analysis and a posteriori error control for the DPG$^*$ method are considered in detail. Reports of several numerical experiments are provided which demonstrate the essential features of the new method. A notable difference between the results from the DPG$^*$ and DPG analyses is that the convergence rates of the former are limited by the regularity of an extraneous Lagrange multiplier variable.

Leszek Demkowicz, Jay Gopalakrishnan, Sriram Nagaraj et al.

A spacetime Discontinuous Petrov Galerkin (DPG) method for the linear time-dependent Schrodinger equation is proposed. The spacetime approach is particularly attractive for capturing irregular solutions. Motivated by the fact that some irregular Schrodinger solutions cannot be solutions of certain first order reformulations, the proposed spacetime method uses the second order Schrodinger operator. Two variational formulations are proved to be well-posed: a strong formulation (with no relaxation of the original equation) and a weak formulation (also called the ultraweak formulation, that transfers all derivatives onto test functions). The convergence of the DPG method based on the ultraweak formulation is investigated using an interpolation operator. A standalone appendix analyzes the ultraweak formulation for general differential operators. Reports of numerical experiments motivated by pulse propagation in dispersive optical fibers are also included.

Federico Fuentes, Leszek Demkowicz, Aleta Wilder

A discontinuous Petrov-Galerkin (DPG) method is used to solve the time-harmonic equations of linear viscoelasticity. It is based on a "broken" primal variational formulation, which is very similar to the classical primal variational formulation used in Galerkin methods, but has additional "interface" variables at the boundaries of the mesh elements. Both the classical and broken formulations are proved to be well-posed in the infinite-dimensional setting, and the resulting discretization is proved to be stable. A full $hp$-convergence analysis is also included, and the analysis is verified using computational simulations. The method is particularly useful as it carries its own natural arbitrary-$p$ a posteriori error estimator, which is fundamental for solving problems with localized solution features. This proves to be useful when validating calibration models of dynamic mechanical analysis (DMA) experiments. Indeed, different DMA experiments of epoxy and silicone resins were successfully validated to within $5\%$ of the quantity of interest using the numerical method.

Brendan Keith, Philipp Knechtges, Nathan V. Roberts et al.

We explore a vexing benchmark problem for viscoelastic fluid flows with the discontinuous Petrov-Galerkin (DPG) finite element method of Demkowicz and Gopalakrishnan [1,2]. In our analysis, we develop an intrinsic a posteriori error indicator which we use for adaptive mesh generation. The DPG method is useful for the problem we consider because the method is inherently stable---requiring no stabilization of the linearized discretization in order to handle the advective terms in the model. Because stabilization is a pressing issue in these models, this happens to become a very useful property of the method which simplifies our analysis. This built-in stability at all length scales and the a posteriori error indicator additionally allows for the generation of parameter-specific meshes starting from a common coarse initial mesh. A DPG discretization always produces a symmetric positive definite stiffness matrix. This feature allows us to use the most efficient direct solvers for all of our computations. We use the Camellia finite element software package [3,4] for all of our analysis.

Federico Fuentes, Brendan Keith, Leszek Demkowicz et al.

This article presents a general approach akin to domain-decomposition methods to solve a single linear PDE, but where each subdomain of a partitioned domain is associated to a distinct variational formulation coming from a mutually well-posed family of broken variational formulations of the original PDE. It can be exploited to solve challenging problems in a variety of physical scenarios where stability or a particular mode of convergence is desired in a part of the domain. The linear elasticity equations are solved in this work, but the approach can be applied to other equations as well. The broken variational formulations, which are essentially extensions of more standard formulations, are characterized by the presence of mesh-dependent broken test spaces and interface trial variables at the boundaries of the elements of the mesh. This allows necessary information to be naturally transmitted between adjacent subdomains, resulting in coupled variational formulations which are then proved to be globally well-posed. They are solved numerically using the DPG methodology, which is especially crafted to produce stable discretizations of broken formulations. Finally, expected convergence rates are verified in two different and illustrative examples.

Federico Fuentes, Brendan Keith, Leszek Demkowicz et al.

A unified construction of high order shape functions is given for all four classical energy spaces ($H^1$, $H(\mathrm{curl})$, $H(\mathrm{div})$ and $L^2$) and for elements of "all" shapes (segment, quadrilateral, triangle, hexahedron, tetrahedron, triangular prism and pyramid). The discrete spaces spanned by the shape functions satisfy the commuting exact sequence property for each element. The shape functions are conforming, hierarchical and compatible with other neighboring elements across shared boundaries so they may be used in hybrid meshes. Expressions for the shape functions are given in coordinate free format in terms of the relevant affine coordinates of each element shape. The polynomial order is allowed to differ for each separate topological entity (vertex, edge, face or interior) in the mesh, so the shape functions can be used to implement local $p$ adaptive finite element methods. Each topological entity may have its own orientation, and the shape functions can have that orientation embedded by a simple permutation of arguments.

Daniele Boffi, Martin Costabel, Monique Dauge et al.

In this paper we prove the discrete compactness property for a wide class of p-version finite element approximations of non-elliptic variational eigenvalue problems in two and three space dimensions. In a very general framework, we find sufficient conditions for the p-version of a generalized discrete compactness property, which is formulated in the setting of discrete differential forms of any order on a d-dimensional polyhedral domain. One of the main tools for the analysis is a recently introduced smoothed Poincaré lifting operator [M. Costabel and A. McIntosh, On Bogovskii and regularized Poincaré integral operators for de Rham complexes on Lipschitz domains, Math. Z., (2010)]. For forms of order 1 our analysis shows that several widely used families of edge finite elements satisfy the discrete compactness property in p-version and hence provide convergent solutions to the Maxwell eigenvalue problem. In particular, Nédélec elements on triangles and tetrahedra (first and second kind) and on parallelograms and parallelepipeds (first kind) are covered by our theory.