{'fraction_non_alphanumeric': 0.03373625910794761, 'fraction_numerical': 0.007602240660405172, 'mean_word_length': 4.659277797640329, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Extended Group Finite Element Method

Interpolation methods for nonlinear finite element discretizations are commonly used to eliminate the computational costs associated with the repeated assembly of the nonlinear systems. While the group finite element formulation interpolates nonlinear terms onto the finite element approximation space, we propose the use of a separate approximation space that is tailored to the nonlinearity. In many cases, this allows for the exact reformulation of the discrete nonlinear problem into a quadratic problem with algebraic constraints. Furthermore, the substitution of the nonlinear terms often shifts general nonlinear forms into trilinear forms, which can easily be described by third-order tensors. The numerical benefits as well as the advantages in comparison to the original group finite element method are studied using a wide variety of academic benchmark problems.


Introduction
A wide variety of mathematical models use partial differential equations (PDEs) to describe physical processes. A powerful tool for the numerical solution of PDEs is the finite element method, see, for example, [10,5] for theoretical aspects and [20] for practical implementation aspects. While linear problems are straightforward and can be solved directly, nonlinear problems present more difficulty. For example, unlike their linear counterparts, nonlinear problems require iterative methods, such as fixpoint iterations or Newton's method. For finite element discretizations, this means that the computationally expensive assembly of the system is required for each iteration. There are many different approaches to improve the efficiency of the assembly process such as parallelization and optimized algorithms using form compilers, see [20]. However, interpolation methods offer a different approach by circumventing the assembly of the nonlinear system in each iteration.
The group finite element formulation, also known as the product approximation technique, is an interpolation method that eliminates the repeated assembly of the nonlinear system, see [6,13] and references therein. Although there is no direct connection, the group finite element method (GFEM) can be seen as a specific application of the interpolation procedure described in [8]. The GFEM replaces groups of variables with an interpolant on the finite element approximation space, which allows the system to be precomputed. While the standard Galerkin approach (SGA) reassembles the system after each update to the approximate solution, the GFEM assembles the interpolated system once before the iterative method is started, see Figure 1. This approach is still an active and relevant research topic. For example, it has been combined with proper orthogonal decomposition in the context of nonlinear model order reduction in [7]. Additionally, a convergence analysis of the group finite element formulation is presented in [4], which also briefly mentions the use of the GFEM in implicit flux-corrected transport schemes, see, for example, [19]. This paper investigates the benefits of extending the group finite element formulation to allow for general finite element approximation spaces for the nonlinear terms. To the best of our knowledge, previous work has only focused on other approaches: For example, a variety of interpolation techniques, such as least-squares fitting, Lagrange or Hermitian interpolation and cubic spline interpolation, for nonlinear functions are investigated in [8]. Alternatively, the GFEM interpolates nonlinear terms onto the underlying finite element approximation space, see [6,13]. Our generalization, which we call the extended group finite element method (EGFEM), allows for more flexibility and broadens the scope of problems that can be treated. Beyond the decrease in the computational effort by eliminating repeated assembly, no sources known to us discuss the unique structure of the interpolated formulation with respect to multi-linear forms. Therefore, an additional aspect that we study is the underlying tensor structure that arises by introducing new variables in place of the nonlinearities. This structure plays an important role in the efficiency of our method, while also showing great promise for model order reduction techniques.
This paper is organized as follows. The mathematical setting and the original group finite element formulation are described in Section 2. We present our extended group finite element method in Section 3 with a focus on the construction and structure of the general approximations as well as the vectorized evaluation of the nonlinear terms. The numerical realizations of our extended group finite element method based on Picard iterations and Newton's method are compared in Section 4. Finally, Section 5 contains a selection of numerical experiments, ranging from the Burgers' equation, which serves as a benchmark problem for the group finite element method, to the p-Laplace equation, which belongs to a class of problems that the original GFEM struggles with.

Background
We present a model problem, which represents a class of nonlinear PDEs, that we consider in this work in Section 2.1. For simplicity, our considerations are restricted to stationary boundary-value problems, time-dependent and/or more general nonlinear problems can be handled similarly. The model problem is discretized using the finite element method, which leads to a nonlinear system that is solved iteratively. Because each iteration requires the costly assembly of the system, we briefly present the group finite element method in Section 2.2. This method, which is also known as the product approximation, interpolates nonlinear terms onto the finite element approximation space so that assembly only needs to be performed once.

Problem Description.
Let Ω ⊂ R d with d ∈ N be an open and bounded domain with sufficiently smooth boundary Γ = Γ D ∪ Γ N such that Γ D ∩ Γ N = ∅. We consider a nonlinear partial differential equation of the form x ∈ Ω, (2.1a) with mixed boundary conditions where n denotes the unit normal vector along the boundary. We assume that the scalarvalued functions in (2.1) ensure the existence of a weak solution u in the trial space V . A weak solution u ∈ V fulfills the variational problem: for all v ∈V . The major difference between the trial space V and the test spaceV lies in the values on the Dirichlet boundary Γ D . In particular, all functions u ∈ V and v ∈V satisfy u| Γ D ≡ u D and v| Γ D ≡ 0. While the Dirichlet boundary conditions are embedded into the trial space, the Neumann/Robin boundary conditions are naturally enforced in the weak formulation (2.2). For more details on the theoretical aspects of nonlinear partial differential equations, we refer the reader to [22,11] and references therein. The approximation of (2.2) is achieved by restricting the infinite-dimensional trial and test spaces V andV to discrete (finite-dimensional) spaces V h andV h , respectively. This leads to the discrete variational problem, i.e., find for all v h ∈V h . For each of the discrete spaces, we have a set of basis functions {φ i } Nu i=1 and {φ i } Nu i=1 that span their respective spaces. The finite element method constructs the discrete spaces by introducing a partition Ω h of the domain Ω and using local function spaces that are defined by a set of finite elements [10,20]. The discrete solutions u h are represented as a linear combination of the basis functions, i.e., Testing (2.3) with each of the test space basis functions delivers the discrete nonlinear problem: for all 1 ≤ i ≤ N u . These nonlinear equations are solved numerically using some form of iterative method, which may either be based on a Picard/fix-point or Newton iteration, see [20] for more details. However, the major difficulty for solving these nonlinear finite element problems lies in the repeated computation of the nonlinear terms, e.g., Because these terms often require numerical integration in order to evaluate, alternative methods that eliminate the repeated assembly of the discrete system are of great importance.

2.2.
Group Finite Element Method. The group finite element method simplifies the nonlinear problems in order to avoid the assembly of the nonlinear system in each iteration of the numerical solution method, see [6,13]. Assume that the given basis {φ i } Nu i=1 is a nodal basis, i.e., φ i (x j ) = δ ij at the degrees of freedom x j associated to the finite element ansatz. This allows groups (products) of variables to be easily interpolated onto the trial space V h , e.g., the scalar-valued nonlinear function f is approximated through which follows from the fact that f (u h (x i )) = f (u i ). The advantage of this method in comparison to the standard Galerkin approach (SGA) is that the nonlinear terms become (multi-)linear forms that can be precomputed, e.g. , Because the integration only needs to be performed once, the computational overhead of each iteration is reduced noticeably. Additionally, the nonlinearity is embedded into the point-wise evaluation of f at the nodes.
Although the GFEM is approximative in nature, there are select examples that display greater accuracy in the nodal values, when compared to the standard approach, see [6,13]. However, these results depend directly on the norm used to compute the error. For one of the examples, the L 2 -norm converges up to two orders of magnitude slower than the standard Galerkin approach, as shown in [6]. Our goal is to extend the GFEM in such a way that the reformulated problem is as exact as possible, while maintaining the reduced complexity and performance gains of the original method. Additionally, we generalize the form of the problems considered in [6] to include the more general form (2.1), which leads to trilinear forms.

Extended Group Finite Element Method
In this section, we extend the group finite element method in order to allow for more generality in the nonlinear terms and their approximations. The latter is achieved by introducing an approximation space W h for the nonlinearity that is tailored specifically to the problem instead of recycling the trial space V h . We discuss the implications and handle the selection of this general approximation in detail for a scalar-valued nonlinear function f in Section 3.1. Afterwards, we briefly present an efficient approach for the point-wise evaluation of f with the help of interpolation operators in Section 3.2. By considering a general model problem, multi-linear forms result after applying the extended group finite element method. In our case, we attain trilinear forms, which are represented through third-order tensors. The resulting tensor formulation is developed and compared to the standard Galerkin approach in Section 3.3.
3.1. Choice of Ansatz Spaces for Nonlinear Terms. The group finite element method interpolates nonlinear functions onto the trial space V h . While this simple choice may be viable in many cases, we introduce a separate approximation space W h that is spanned by the nodal basis Interpolating the nonlinear function f onto W h instead of V h delivers the following approximation: The coefficients f = (f 1 , . . . , f N f ) T are defined through the point-wise evaluation of the function f on the degrees of freedom, i.e., (3.6) which is motivated by the fact that f h (x k ) = f k . A nice feature of our extended group finite element method is that the original formulation is recovered for W h = V h . However, there are many situations, where a different approximation space W h = V h may be more natural or even necessary in order to avoid an ill-posed problem.
While the coefficients f are defined through (3.6), the approximation space W h still needs to be fixed. In particular, the actual setting, e.g., the nonlinear function f , the trial space V h and other numerical considerations, influences our choice for the local finite element space. The plausible scenarios that influence how to approximate f are separated into the following three cases.
Case 1 (Exact Reformulation). Let T denote the an element in Ω h . At this point, a trial space V h based on finite elements has been selected. When we consider the restriction of f to the element T , it follows that f | T (u) = f (u| T ). This case assumes that the restriction can be described exactly by some finite element. This leads us directly to an appropriate approximation space W h .
Example 3.1. Let P k denote the space of polynomial functions of at most degree k and T an element in Ω h . Consider u h | T ∈ P k . All monomials f (u) := u with ∈ N 0 fulfill f | T ∈ P ·k . Alternatively, the nonlinear function defined through f (u) := ∇u 2 can be characterized locally through f | T ∈ P 2 (k−1) . Case 2 (Approximation). Let the trial space V h and the approximation space W h be fixed. Then, the nonlinear function f is approximated locally by the finite elements used in the construction of W h . This case depends on a variety of factors. For example, if Case 1 is not applicable, then it may be acceptable to approximate the nonlinearity directly. Additionally, some problems require other approximations, e.g., errors that arise from evaluating the integration in (2.3) inexactly. The current case allows for a direct approximation of f through a given W h , which takes these situations into consideration.
Example 3.2. The group finite element method falls into this case. Here, the approximation space is the same as the trial space, i.e., V h = W h . Case 3 (Quadrature). For a fixed trial space V h and a single element T in Ω h , we consider the situation, where an exact reformulation of the discrete problem is desired but not feasible with existing finite element spaces. In particular, we construct an approximation space W h by defining a set of discontinuous basis functions {η k } N f k=1 that span it. For this purpose, we propose a "new" type of discontinuous finite element, which is based on the quadrature rule used to numerically integrate over the element T . For example, a given q-point quadrature rule approximates the integral where w denotes the weights of the quadrature rule and x denotes the quadrature nodes located on the element T . In order to define a finite element, we fix the degrees of freedom and basis functions on an element. Here, the degrees of freedom are given by the quadrature nodes {x } q =1 . The basis functions are defined through indicator/Dirac-delta functions related to the quadrature nodes, i.e., Using these basis functions, the nonlinear function f is approximated by where f is defined through (3.6). Introducing this approximation into the above integral leads to . Therefore, the quadrature is embedded into the approximation space W h and the discrete version of the original problem is conserved. Note that W h is constructed by combining the quadrature points {x } q =1 and basis functions {η } q =1 for each element T . Depending on the chosen quadrature rule, this quickly leads to a large number of unknowns. For example, the quadrature points of symmetrical Gaussian quadrature rules on triangles are located in the interior of the triangle [9]. Therefore, the approximation space constructed using such a quadrature rule consists of q N el basis functions and degrees of freedom, where q and N el denote the number of quadrature points and the number of elements in the mesh Ω h , respectively. Example 3.3. Consider a triangular element T k from Ω h . The following quadrature rule is exact for linear functions: where c k denotes the center point of the triangle and |T k | the triangle's area. The associated quadrature finite element consists of a single degree of freedom at the center c k with the constructed basis function η k defined in (3.7) with w k = |T k |. Compared to the discontinuous constant finite element, the location of the degree of freedom coincides, while the basis functions are completely different.

Algebraic Equations.
After introducing the approximation space W h for the nonlinear function f in Section 3.1, the coefficients f are defined through the point-wise evaluation of the function f . These algebraic equations require the interpolated values and/or gradients of the approximate solution u h at every degree of freedom belonging to f , i.e.
. Instead of evaluating each constraint individually, the evaluation can be vectorized with the help of two interpolation operators. First, the evaluation of the approximate solution u h on each degree of freedom x i associated to f can be performed at once using Similarly, a third-order tensor interpolates the gradient of the approximate solution ∇u h onto the degrees of freedom x : This leads to the "vectorization" of the algebraic equations via , which allows for the efficient evaluation of these terms in the numerical implementation.

3.3.
Multi-linear Forms. The idea of introducing variables as placeholders for nonlinear terms in order to introduce structure to a given problem can be found in many different applications. For example, the idea of quadratization [15,18], which reformulates nonlinear problems as quadratic-bilinear systems, can be found in the field of model order reduction. Alternatively, "lifting" approaches for Newton's method bring structure to nonlinear optimization problems, which often converge faster than the nonlifted variants [1]. Similarly, the group finite element formulation and our extended variant also introduce a clear structure to the discrete nonlinear problem (2.4). Here, the nonlinear terms in (2.5) simplify into multi-linear forms. Note that depending on the nonlinear functions, different forms may arise but can be handled similarly. Our focus in this subsection lies in the structure introduced by interpolating the nonlinearities onto approximation spaces. In particular, nonlinear forms are replaced with multi-linear forms. These simpler forms are represented by appropriate tensors, which prove to be numerically much more efficient in comparison to the assembly of the nonlinear terms.
Notation 3.4. Let T denote a third-order tensor, A a second-order tensor (matrix) and v, w first-order tensors (vectors). We maintain this notation for the remainder of this work. Then, the tensor (dyadic) product is defined through The first and second contractions, denoted by · and :, are defined through Additionally, the first contraction can also be defined along a different dimension, which we denote by appending an index to the ·, e.g., as already introduced in Section 3.2. Finally, we denote the Hadamard product, which refers to an element-wise product, through for vectors of the same dimension. From these definitions, it is clear that the tensor product results in a second-order tensor, while the contractions reduce the order of the tensor by one and two orders, respectively. The Hadamard product does not change the order/dimension at all.
In order to facilitate the discussion regarding the tensor structure of the extended group finite element formulation, we focus on the discrete problem (2.4). For the standard Galerkin approach, which requires the assembly of the nonlinear terms in every iteration of the iterative method, we define the following notation using the relationship where the forms associated to K and M are similar to the classical stiffness and mass matrices. The boundary integral in (2.5) is handled similarly to the second term, requiring only a different integration domain, and is denoted by g. In contrast to the standard formulation, our extended group finite element formulation has a well-defined structure that can be precomputed. Applying the discretization discussed in Section 3.1 to (2.5) results in the following discrete forms: The interesting development in comparison to the GFEM is the tensor structure that results from approximating K(a, b, u). By interpolating both nonlinear functions a and b, we attain a trilinear form, which is represented with the help of a third-order tensor.
Remark 3.5. In addition to the tensor notation, it is also possible to express the action of a trilinear form in terms of a matrix-vector product by using matricization and the Kronecker product. In particular, we can rewrite the second contraction with the tensor product as where K b a denotes the mode-2 matricization of the tensor K b a and ⊗ K the Kronecker product. However, due to the size of the systems that result from the mode-2 matricization and Kronecker product, the tensor formulation is much more efficient for the purpose of numerical implementation.

Numerical Methods
The advantages of the extended group finite element formulation become even more clear, when the discrete systems, which must be solved in each iteration, are compared to the standard Galerkin approach. In practice, there are two common iterative methods for solving nonlinear systems such as (2.4), see [20]. The first method uses the Picard iteration and is also known as the method of successive iterations. The idea is to substitute the previous iteration into the nonlinear terms, which delivers a linear system in terms of the new iteration. The second method, Newton's method, is a well-known method for solving nonlinear systems, which exploits gradient information in order to improve the (local) convergence rate. However, this comes at the cost of evaluating the Jacobian matrix, which is also assembled using numerical integration. While Newton's method generally converges more quickly, the Picard iterations are simpler to implement. In light of its simplicity, the benefits of the extended group finite element method are demonstrated using the Picard iterations in Section 5. The results are qualitative and, therefore, expected to be similar when using Newton's method to solve the nonlinear systems.
The standard Galerkin approach solves the following nonlinear form with the nonlinear terms defined in (3.8), while the extended group finite element formulation leads toF with z :=(u, a, b, c, g) T and the new interpolated terms defined in (3.9). The point-wise definitions of the coefficient vectors a, b, c and g assume that the associated function evaluations of a, b, c and g are vectorized, e.g., which is realized with the help of the interpolation operators Π f u and Π f ∇u from Section 3.2, for f ∈ {a, b, c, f }.
Picard Iteration. The form of the Picard iterations is problem-dependent. In general, the desired form of the discrete nonlinear problem is given by In this setting, the Picard iteration is attained by substituting the previous iteration into the nonlinear terms, i.e., u n+1 = A(u n ) −1 r(u n ).
In order to illustrate this idea for a simplified version of (2.1), let b(u) = u. Then, the original term K(a, b, u) in (4.10) can be written as K(a, u) u with (K(a, u)) ij := Ω h a(u h (x)) ∇φ j (x) · ∇φ i (x) dx, while the tensor form K b a : (b ⊗ a) in (4.11) simplifies to K a : (u ⊗ a) = (K a · a) · u with This leads to the following iteration for the standard Galerkin approach: In contrast, the (extended) group finite element formulation iterates over Although these forms are only valid for the case that b ≡ id, they still illustrate that the reformulated systems avoid assembly in each iteration. Instead, the nonlinear functions are evaluated point-wise, which is computationally less expensive than the numerical integration of the nonlinear forms.
Newton Iteration. The iterations defined by Newton's method can be attained through two different approaches [20]. The standard approach is to apply Newton's method to the discrete nonlinear problem (2.4), while alternatively it can also be applied to the continuous nonlinear PDE (2.1). In the end, both approaches lead to the same system. The iterates are given using the following update: where D u F denotes the Jacobian matrix. For the standard Galerkin formulation, this system uses (4.10) with the Jacobian matrix .
In order to evaluate the Jacobian matrix, assembly using numerical integration is necessary. This leads to a noticeable increase in the computational complexity of each iteration in comparison to the Picard iterations although the Newton iterations are expected to converge more quickly. The extended group finite element formulation iterates over z = (u, a, b, c, g) T using withF defined in (4.11) and the Jacobian matrix In contrast to the standard Galerkin formulation, the Jacobian matrix here does not require any additional integration. Instead, the only integrated portions depend on the precomputed forms K b a , M c and G g . In addition, the Jacobian matrices D u a, D u b, D u c and D u g associated to the nonlinear functions a, b, c and g are defined point-wise.
Remark 4.1. In some cases, it may be necessary to use some kind of dampening or continuation method in order to ensure convergence of the iterations. Because this depends entirely on the nonlinearities and the formulation of the discrete problem, there is no general guideline to determine if such methods are necessary. An example of such a problem is presented in Section 5.3.

Results
In this section, we present a comprehensive performance study of the extended group finite element method (EGFEM) for a variety of examples. The examples serve to investigate different aspects of the extensions presented in Section 3, especially in comparison to the standard Galerkin approach (SGA) and original group finite element method (GFEM). In addition to the advantages and disadvantages of the tensor formulation, the exact reformulation of the EGFEM contrasts the approximative nature of the GFEM. Each example begins with a brief description of the nonlinear problem, where the structure is presented with respect to the general form (2.1). Afterwards, the discrete systems related to the SGA as well as the original and extended GFEMs are identified. Finally, the important aspects of the given example are discussed and supported with the help of figures and tables. The wide variety of the examples shows the versatility and applicability of our extended group finite element method for nonlinear problems.
All simulations are performed in Matlab * using Sandia's Tensor Toolbox for the tensor calculations and Gmsh for the mesh generation [3,2,14]. For simplicity, the solution u is * All results were computed on a machine with an i7-8700 and 32 GB of RAM using Matlab R2019b. discretized using linear finite elements in all of the examples. The assembly of the system is performed with the help of symmetrical Gaussian quadrature rules for triangles [9]. In light of the additional complexity associated to Newton's method, we use the Picard iterations in order to solve the nonlinear problems to a given tolerance. Also note that the runtimes presented refer only to the online computational times, i.e., the time spent in the iteration method. The time spent computing the system tensors, matrices and vectors as well as the interpolation operators is neglected, since most relevant applications such as many-query and/or real-time optimal control perform these computations in an offline phase. However, it should be noted that these calculations can be time-consuming. 5.1. Quadratic Nonlinearity. A simple example of a nonlinear problem is the Poisson problem with an additional quadratic term: The source term d is chosen according to the method of manufactured solutions [21], where the desired solution is given through the extension of u D onto the entire domainΩ, i.e., In relation to (2.1), this example is given through with Γ D = Γ and Γ N = ∅. After discretizing with linear Lagrange elements, the standard Galerkin approach delivers the following nonlinear discrete problem: with the stiffness matrix K. An interesting aspect of this example is that the nonlinearity c can be handled in two different ways. Firstly, it can be treated directly as the trilinear form, e.g., This eliminates the need for repetitive assembly of the nonlinear term in each iteration, i.e., the discrete problem is given through Alternatively, the (extended) group finite element method introduces the variable c, which is defined point-wise. This leads to the replacement of the second term M(c, u) and introduction of an additional equation , where the interpolation operator Π c u is equal to the identity in the case of the original group finite element formulation. While the GFEM chooses W h = V h , our EGFEM allows for a more general choice. The space of piece-wise quadratic functions (P 2 ) exactly describes the nonlinearity. Hence, the nonlinear function c falls into Case 1 of Section 3.1 for the extended group finite element method with quadratic Lagrange elements for W h . For the sake of comparison, we also consider Case 3 by using a quadrature rule that is exact for cubic polynomials (I 3 ).
A comprehensive overview of the numerical performance of each approach can be found in Figure 2. When comparing the runtime and, therefore, also the speedup, we see that the tensor formulation becomes more efficient as the size of the problem grows. However, the original and extended group finite element formulations immediately show an improvement  Figure 2. An overview of the performance of each approach for the quadratic example using a regular mesh that is successively refined through splitting. From left to right: The amount of time to solve the system for increasingly finer meshes, the speedup describes the ratio between the reference runtime (SGA) and the considered approach, and the relative error with respect to the approximate L 2 -norm from the exact solution.
that delivers a speedup between one to two orders of magnitude. The important difference between the GFEM and the EGFEM can be seen in the relative error with respect to the exact solution. Here, the GFEM converges slower and with a larger error than the extended formulation. This lies on the approximative nature of the group finite element formulation. In particular, the GFEM interpolates the nonlinearity onto the space of piecewise linear functions (P 1 ). By choosing the space of piece-wise quadratic functions (P 2 ), the EGFEM preserves the nature of the nonlinearity. Therefore, the extended group finite element formulation is numerically equivalent to the standard Galerkin approach, while exhibiting the same increase in performance as the original formulation.

Burgers' Equation.
A common example used in many literature sources for illustrating the group finite element formulation is Burgers' equation, see, for example, [6,13,7]. Here, we consider the viscous two-dimensional Burgers' equation with homogeneous Dirichlet boundary values on the time interval [0, T ] with T > 0, i.e., with the initial value u(x, 0) = u 0 (x) for all x ∈ Ω and u(x, t) = 0 on the boundary. Adopting the same setting as [7], the initial value u 0 and the source d are chosen so that the exact solution is given through the following function: while we choose ν = 1 and T = 1. By taking advantage of the fact that the nonlinear terms can be rewritten as we apply integration by parts in order to attain the following semi-discrete form for the standard Galerkin approach: However, since the action N(u(t)) can also be described with the help of a multi-linear form, we also introduce the following tensor-based formulation, which does not require repeated assembly of the system: where the tensor N is defined through The extended and original group finite element methods interpolate the term u 2 onto the approximation space W h and V h , respectively. This leads to the semi-discrete form We apply a semi-implicit discretization for the temporal component in order to attain an iteration rule similar to the Picard iterations, i.e., by using the previous iterate for the nonlinear terms and the new iterate for all other terms with the time step δt. This leads to for the standard Galerkin approach. By using a tensor representation for the multi-linear form, the SGA can also be defined through The (extended) group finite element formulation solves the following system in each time step:F The equation for u n+1 is linear, since the nonlinear terms are evaluated using the previous iteration. Afterwards, the nonlinear coefficient f n+1 can be updated with the help of the new iterate u n+1 . Figure 3 compares the performance of each method, when using δt = 10 −2 for the time discretization. We see that the tensor-based formulation improves the performance, since the computational overhead associated to assembling the system in each iteration is removed. Additional performance gains are achieved by introducing the (extended) group finite element formulation, which relies on matrices instead of a third-order tensor for the nonlinear terms. This difference leads to the slight improvements between the (extended) group finite element approach and the tensor formulation. In contrast to the first benchmark problem, we see that the difference in accuracy between the original and extended group finite element formulations is numerically negligible. While the GFEM interpolates the nonlinearity onto a piece-wise linear approximation space (P 1 ), the EGFEM exactly describes the nonlinearity in the discrete setting by using the finite element space composed of piece-wise quadratic functions (P 2 ). When compared to the exact solution, the error differences are minimal. Based on these and following results, we see that Case 2 in Section 3.1 is suitable in many different situations. However, as we have also seen in Section 5.1, the EGFEM can lead to noticeably more accurate results with similar  Figure 3. An overview of the performance of each approach for the viscous Burgers' equation using a regular mesh that is successively refined through splitting. From left to right: The average amount of time to solve the problem for increasingly finer meshes, the speedup describes the ratio between the reference runtime (SGA) and the considered approach, and the relative error with respect to the approximate L 2 -norm from the exact solution.
performance gains. In the end, the choice of the approximation space depends on the setting and problem to be solved.

Superconductivity.
For this example, we consider a simplified Ginzburg-Landau model for superconductivity as presented in [16,17]. The semi-linear model is given by x ∈ Γ.
We use the method of manufactured solutions in order to determine the source d as The Dirichlet boundary condition is chosen as u D (x) = 1 6 sin(2 π x 1 ) sin(2 π x 2 ) exp(2 x 1 ).
In the framework of (2.1), we have a(x, u, ∇u) = ν, b(u) = u, c(x, u, ∇u) = u 3 + u with Γ D = Γ and Γ N = ∅. This example allows us to investigate two different aspects, which can influence the convergence of the iterations. On the one side, the parameter ν controls the importance of the nonlinearity c for 0 < ν 1. On the other hand, the nonlinearity c can be decomposed in the following three different ways: where the (weighted) mass matrix M(c, u) is defined through In total, we consider the following three systems for the standard Galerkin approach, which are all analytically equivalent: with the respective nonlinearityc. Each way of decomposing the nonlinear function c results in different discrete systems with different unknowns. This influences the choice of the approximation space W h , when applying the extended group finite element formulation. We consider a range of choices, including piece-wise quadratic functions (P 2 ) and a quadrature rule that is exact for quartic polynomials (I 4 ) for (5.13a) and piece-wise cubic functions (P 3 ) as well as the same quadrature rule as before for (5.13b) and (5.13c). The group finite element method chooses W h = V h in each situation. The discrete problems solved by the (extended) group finite element method are given through with appropriate point-wise definitions for the auxilliary variable c orc. For the numerical tests, we consider ν ∈ {1, 10 −2 , 10 −3 }, since for decreasing values of ν the nonlinear term becomes critical. A complete overview of these results for a fixed mesh can be seen in Table 1. The empty entries of the table denote the formulations, where the fixed point iteration either diverges or oscillates. For the situations in which the method converges, we see a noticeable decrease in the computational effort for the reformulated problems. Figure 4 depicts the performance of the discussed methods for ν = 1 and the setting (5.13a). Here, we see the improvements in performance for the systems that avoid repeated assembly of the nonlinear forms. The error in Figure 4 is dominated by the approximation error, which leads to only a slight variation between the reformulated systems. The deviation from the error associated to the standard Galerkin approach shows the error introduced by approximating the nonlinearity with the GFEM in comparison to the exact reformulation of the EGFEM. In summary, a careful balance between the size of the reformulated system and the approximation errors must be found. In situations such as this example, Case 2 in Section 3.1 -only approximating the nonlinearity -are viable and more efficient than the exact alternatives.
A very interesting aspect of the reformulation process becomes apparent in Figure 4. From the absolute errors and the similarity between each method, we infer that the approximation error from u h dominates the total error. This means that the approximation error from the nonlinearity becomes negligible, as seen in the deviations shown in Figure 4. In light of this, the second approach for handling the algebraic constraints in Section 3.1 may be viable in many situations, where this phenomenon occurs. In this case, the increased computational effort for exactly reformulating the nonlinear terms is unnecessary and can be avoided.    Figure 4. An overview of the performance of each approach using (5.13a) with ν = 1. The source term is assembled using an adaptive quadrature rule, which reduces the numerical errors from the numerical integration as much as possible. From left to right: Runtime in seconds to solve the problem, speedup associated to each formulation (with respect to the standard approach) and the relative error with respect to the exact solution using the approximate L 2 -norm.
Similar to Section 5.3, we have a situation, where the nonlinear term c can either be written as However, in light of the previous results, we only consider the first case. Here, we have a situation, where the nonlinearity -a rational function in terms of u -does not have a Abs. Error Figure 5. Convergence of quadrature rule. From left to right: absolute error between the solutions using first and second order, second and third order and third and fourth order quadrature rules, respectively. matching finite element. Therefore, we consider Case 2 and Case 3 of Section 3.1. Using the (extended) group finite element formulation, the system can be rewritten as The only problem that remains is which quadrature rule makes the most sense in order to evaluate the nonlinear term. If we look at the original term and use the information that u h | T ∈ P 1 , then we can see that where P , Q are at most polynomials of degree 2 and 1, respectively. Using this knowledge, it makes sense that the rational term P/Q behaves at most quadratically, and, therefore, a quadrature rule, which is exact for quadratic polynomials, delivers sufficient results. This hypothesis is confirmed in Figure 5, where quadrature rules of increasing order are compared.
The results shown in Figure 6 lead to interesting conclusions. First, we see that the quadrature elements ensure that the standard approach is conserved. In contrast, the approximations, while not being exact reformulations of the original problem, offer excellent results even with piece-wise constant elements. However, the absolute error is dominated by the approximation error resulting from the discretization from u h and not the approximation of the nonlinearity. The deviations in Figure 6 show the influence of the reformulation technique. Interestingly, the higher order approximation using quadratic finite elements leads to a more accurate approximation in the coarser meshes. This discrepancy can also be associated with the higher order quadrature rule, which is required to assemble  Figure 7. An overview of the performance of each approach for the p-Laplace equation using a regular mesh that is successively refined through splitting. From left to right: The average amount of time to solve the problem for increasingly finer meshes, the speedup describes the ratio between the reference runtime (SGA) and the considered approach and the relative error with respect to the L 2 -norm from the exact solution.
the resulting tensor, see Figure 5 for the influence of the quadrature rules in the solution accuracy.

p-Laplace Equation.
The p-Laplacian, which generalizes the Laplacian operator, defines a quasilinear elliptic partial differential operator. The p-Laplace equation with Dirichlet boundary conditions is given through For the case d ≡ 1, u D ≡ 0 and Ω = {x ∈ R 2 : x ≤ 1}, an explicit solution is given by for p ∈ (1, ∞). In our setting of (2.1), we have a(u) = ∇u p−2 , b(u) = u, c(u) = 0, d = 1 with Γ D = Γ and Γ N = ∅. The discretization using linear finite elements leads to the Picard iterations for the standard approach K(a, u n ) u n+1 = d and for the extended group finite element formulation K : (u n+1 ⊗ a n ) = d, a n+1 = a(Π a ∇u · 2 u n+1 ). Because the gradient ∇u h is discontinuous along the interfaces of the elements, it is not possible to interpolate the nonlinearity a onto the trial space V h . Hence, the GFEM is not applicable. For this example, we consider the EGFEM with W h using piece-wise constant elements (P 0 ) and the midpoint quadrature rule (I 1 ), which is exact for linear functions.
The computational advantages of the extended group finite element formulation are depicted in Figure 7. Additionally, the extended formulations are equivalent with the standard Galerkin approach. In such situations, performance gains are achieved at the cost of increasing storage requirements, e.g., for storing the third-order tensors, which can be very large but are also sparse, and the intermediate variables, which store the information about the nonlinearity.  Figure 8. An overview of the performance of each approach for the minimal surface equation using a regular mesh that is successively refined through splitting. From left to right: The average amount of time to solve the problem for increasingly finer meshes, the speedup describes the ratio between the reference runtime (SGA) and the considered approach, and the relative error with respect to the approximate L 2 -norm from the exact solution.

Minimal Surface Equation.
For this example, we consider the minimal surface equation, which is a quasi-linear partial differential equation, i.e., If the source d vanishes, then the solution describes the function whose graph has the minimal area of all graphs over Ω with the fixed values on the boundary [12]. We choose the unit square [0, 1] × [0, 1] as the domain and define the boundary condition as for all x ∈ Γ. For the numerical investigations, the source term d is chosen according to the method of manufactured solutions [21], where the desired solution is given through the extension of u D onto the entire domainΩ. This allows us to study the accuracy of the method by comparing the numerical solutions to the exact solution, while also comparing the performance of the extended group finite element formulation with the standard Galerkin approach. In relation to (2.1), this example is given through a(x, u, ∇u) = (1 + ∇u 2 ) − 1 2 , b(u) = u, c ≡ 0, d = −∇ · ∇u D (x) 1 + ∇u D (x) 2 with Γ D = Γ and Γ N = ∅. After discretizing with linear Lagrange elements, the standard Galerkin approach delivers the following nonlinear discrete problem: K(a, u n ) u n+1 = d.
In contrast, the extended group finite element method solves K : (u n+1 ⊗ a n ) = d, a n+1 = a(Π a ∇u · 2 u n+1 ). In this situation, the GFEM is not viable, since the gradient of the linear approximation u h is not defined on the nodes. Therefore, the extension presented in this work is necessary. In light of the fact that a is constant on each element, we consider the cases that W h is given through piece-wise constant functions (P 0 ) and using a quadrature rule that is exact for linear functions (I 1 ). Figure 8 illustrates the advantages of the EGFEM. The average runtime and, consequently, the respective speedup testify to the performance increases by eliminating the repeated assembly of the nonlinear problem. For this benchmark problem, the reformulation using intermediate variables for the nonlinearity is numerically exact as can be seen in convergence of the absolute error of the solutions under mesh refinement. In summary, the extension to more general approximation spaces W h is necessary in order to apply the group finite element formulation to this problem, which results in a numerically equivalent structured problem that can be solved faster than the standard approach.

Discussion
The assembly of the nonlinear systems accounts for a large amount of computational overhead. Therefore, techniques like the group finite element formulation are crucial for time-sensitive computations. An important and relevant example is the optimal control of a nonlinear PDE system. Because gradient-based methods require the repeated solution of the forward and backward problems, the original and our extended group finite element methods become invaluable, since large portions of the systems only need to be computed once. However, the original group finite element formulation has its limitations. By interpolating all nonlinear terms onto the trial space V h , the resulting system is only an approximation of the original problem. At the same time, many difficulties arise for nonlinear problems that depend on gradient information. In such situations, the original group finite element formulation may not even be applicable. In order to allow for more flexibility as well as improve the accuracy of the reformulated systems, we extend the group finite element method to use general finite element approximation spaces. As seen in Section 5, this allows for the exact reformulation of the original problem in many cases as well as the handling of nonlinear functions that depend on the gradient of the approximate solution u h . The additional extension of the group finite element formulation to multi-linear forms using tensors introduces structure to problems that can be used even without the introduction of intermediate variables, see, for example, Sections 5.1 and 5.2. These adaptations lead to a more efficient approximation of nonlinear problems with the same computational benefits as the group finite element method. While our results use the Picard iterations, it is reasonable to expect qualitatively similar results when using Newton's method. This is because the most important aspect of the extended group finite element method is the elimination of the computational overhead associated with repeated assembly. A real-world application using our proposed method is presented in [23]. Additionally, further extensions including the combination with model order reduction techniques, similar to [7], are currently being investigated.