Error estimation and adaptivity for PGD based on complementary solutions applied to a simple 1D problem

Reduced order methods are powerful tools for the design and analysis of sophisticated systems, reducing computational costs and speeding up the development process. Among these reduced order methods, the Proper Generalized Decomposition is a well-established one, commonly used to deal with multi-dimensional problems that often suffer from the curse of dimensionality. Although the PGD method has been around for some time now, it still lacks mechanisms to assess the quality of the solutions obtained. This paper explores the dual error analysis in the scope of the PGD, using complementary solutions to compute error bounds and drive an adaptivity process, applied to a simple 1D problem. The energy of the error obtained from the dual analysis is used to determine the quality of the PGD approximations. We define a new adaptivity indicator based on the energy of the error and use it to drive parametric h- and p- adaptivity processes. The results are positive, with the indicator accurately capturing the parameter that will lead to lowest errors.


Introduction
Computational methods have become an essential part of most high end engineering projects.They greatly simplify the design and analysis of highly complex systems and are a must for companies that aim to become and stay competitive.When modeling most demanding systems a high computational effort is needed, leading to a low response time between the identification of the details of the model to be considered and the availability of the response.Furthermore, the modification of details in the model leads to re-analyses of the process, further delaying the determination of the desired response.This is a major issue in our current quick paced reality, which relies in the premise that the faster you can come up with a reliable solution, the better.Reduced order methods appear as an answer for this demand.
The main idea behind reduced order methods is to formulate a model that takes only the essential parts of a simulation, reducing the computational time needed to perform Fig. 1 Linear elastic straight bar on an elastic support, with n b sections.The bar is fixed at the start and has a concentrated load at the end a complex analysis while aiming to maintain the accuracy of the results.The Proper Generalized Decomposition (PGD), which we use in this work, is one of the several reduced order methods that exist.The main characteristics of this method is its a priory nature, which avoids the need to perform the full simulation beforehand in arbitrarily selected instances [1].
One of the challenges that the application of PGD faces is that the method lacks a posteriori estimations tools and adaptivity strategies [2][3][4][5].Another issue with PGD is its complexity when dealing with geometry as one of the parameters [6,7].
This paper has the objective of showing the application of a PGD driven adaptivity process to a simple 1D problem.The particular aspect of this implementation is that we simultaneously seek two complementary PGD approximations, one compatible and one equilibrated, which we will use to bound their error (as in [8]) and also to drive the adaptivity process (in the physical and in the parameter space).
The governing equations that describe the problem being considered are presented for the case of a bar divided in two sections, followed by a discussion on error assessment strategies to evaluate the solutions obtained.We proceed to explain the steps used to obtain approximated solutions, first for the finite element method and later using the proper generalized decomposition.We present the parametric form for the approximated solutions and a form of assessing the PGD solutions error, this time specific for each parameter.We apply this error measurement strategy to obtain a novel error indicator which can drive both h-and p-adaptivity processes.

Governing differential equations
We consider a simple problem for which an analytic solution can be obtained: a linear elastic straight bar divided in n b sections, each section with a length (γ ), an uniform axial stiffness (EA) and an uniform elastic support (k), connected in sequence, which form a one dimensional lattice structure.Note that when the stiffness of the support is zero, the solution of this problem is trivial.
Although arbitrary imposed displacements or applied forces may be considered at the nodes that connect the different sections, we restrict our study to the problem where the displacement ( ) is fixed at the start of the first bar and a non-zero force (P) is considered at the free end of the bar, as presented in Fig. 1.
We want to obtain a solution for this problem, expressed in terms of either the displacements (u) or of the axial forces (N ), that is characterized by the values of γ , EA and k for each section, plus the value of the imposed displacement and of the applied force.This results in, at most, (3 × n b ) + 2 parameters.
We know that the strain in an arbitrary point at the bar is The corresponding axial force is where E is the Young's modulus at the specified point and A its cross section area.Assuming that the force transmitted by the support is F = ku, with k being the elastic stiffness of the support, and knowing that this force has to be balanced by the axial force, we have: The exact solution of the problem has to satisfy: subjected to EA du dx x=L = P. ( Notice that for a bar divided in multiple sections we will need to impose the continuity of displacements and of the corresponding axial force in the n b − 1 connections nodes of the sections.We can also write the solution for this problem in terms of the axial force, considering first equilibrium and then compatibility.We know that for a given axial force distribution, equilibrium will be satisfied if relation (3) is true, implying that: Compatibility and the constitutive relations require that the strain corresponding to this displacement has to be equal to the strain associated with the axial force, such that: If we write the problem in terms of the axial force, which has to be continuous, it becomes subjected to Again, we need to impose the continuity of the axial forces and of the corresponding displacements at the connection nodes of a bar divided in multiple sections.

Error assessment measurements
In this section we discuss two possible ways to compute the error for approximate solutions.This is made for an arbitrary problem with domain and boundary , which for the bar that we are considering are is its length and the endpoints, respectively.We will first talk about the error in energy and later of the energy of the error.After a brief explanation, these errors indicators will be applied to our problem.Both the error in energy and the energy of the error can work as global estimators of the solution's behavior, providing error bounds.

Error in Energy
Following [9], the potential energy of a mechanical system can be defined as a function of a kinematically admissible displacement u, such that: with U describing the total strain energy, W the strain energy density, and V is the work done by the applied forces, defined as: where b are the body forces, t the boundary traction and t the static boundary.The complementary potential energy c for a statically admissible stress field N can be written in a similar way as where U c is the complementary strain energy, W c the complementary strain energy density, and V c is the work done by the imposed displacements.We can define V c as with u representing the kinematic boundary, N the boundary operator, which for the 1D problem is ±1, and ũ the displacement at the boundary.Considering linear elastic constitutive relations and that the possible influence of initial strains is excluded, we can write the strain energy density as: and the complementary strain energy density: We can use the dual solutions to write the average values for the total potential energies and the strain energies, such that: The sum of the potential energy and the complementary potential energy is zero for the exact solution [9].For force driven or displacement driven problems, such as the one that we will consider, the error of each approximate solution is bounded by the sum of their energies.Compatible solutions for force-driven problems have the total energy equal to minus the strain energy, which is a negative value.We choose to work with the symmetric of the total energy, which is positive and converges from bellow, and can be used as an error measure.On the other hand, compatible solutions for displacement-driven problems will present a total energy which is always equal to the strain energy, therefore converging from above.Equilibrium solutions will have an analogous opposite behavior, where forcedriven problems will converge from above and displacement-driven from bellow (More about this in [8]).

Energy of the error
We can compute the exact energy of the error of a compatible or of an equilibrium solution, respectively as: where the subscript k indicates solutions that come from kinematically admissible displacements, s solutions from statically admissible axial forces and the lack of a subscript indicates the exact solution.An upper bound of both values is [10]: Unlike the energies of the solutions, this expression is applicable to any type of problem (force driven, displacement driven or mixed).Therefore, from now on we always use the energy of the error instead of the error in an energy.We call the integrand in (19) an error density, defined as: Therefore, we can rewrite Eq. (19) as: It can be shown that [9]:

Compatible and equilibrated finite elements approximations
For our problem, any general function f (x) ∈ H 1 (continuous and differentiable) that satisfies f (0) = is a compatible displacement field, with a normal force distribution N = u EA.Also, any general function g(x) ∈ H 1 is an equilibrated normal force distribution, provided g(L) = P, with a displacement field u = g /k.
When considering the displacement field, we can write the bilinear form of the strain energy product as: and the work of the concentrated force is where a k (u α , u β ) and b k (u β ) are the bilinear and linear form products, corresponding to compatible (kinematically admissible) solutions.
For the equilibrated solution (statically admissible), we will have: and the work of the imposed displacement is It is also possible to consider a mixed form We use a Galerkin approach to determine the approximate solutions.Assuming a polynomial degree of approximation p, we have: where vectors û and N represent the coefficients of the approximations for the displacements and axial forces, respectively, and φ(x) are the polynomial basis.Substituting (28) into ( 23) and (25) will lead to the definition of the stiffness and flexibility matrices, K and F , respectively, (29) We can define a matrix that projects the equilibrated solutions onto the compatible ones, as in (27), such that: The vectors of applied forces and imposed displacements are obtained substituting (28) into ( 24) and (26), so that: The resulting finite element compatible and equilibrated systems are, respectively We can obtain the finite element approximation for the displacements u and axial forces N by solving the equations in (33) and using the results at Eqs. (28).

Compatible and equilibrated PGD approximations
A PGD formulation is going to be used to obtain the approximations of the solutions of each complementary problem, as in [8].We define μ, a vector of D parameters μ 1 , μ 2 , . . . ,μ D , each defined in its own domain μ i ⊂ R, with i = 1, 2, . . . ,D. Then, the vector of parameters μ is defined so that μ Each solution will depend on the D parameters and will be approximated as a sum of N modes of D independent functions of each variable.We can write a general PGD approximated solution as: with f representing the coefficients of the approximation in the physical space.A linear combination of polynomials of degree p j is used to represent the function: which is a generalization of (28) where F represents the coefficients of the approximation in the parameters.
A fixed point iteration is used to determine the terms in each product, or more precisely, the coefficients F and f .This iteration procedure will try to minimize the extended Galerkin residual of (23) or (25).For a generic term n > 1, f , is written as we can further simply f as: The extended residual is obtained by considering the integration domain of ( 23) or ( 25) including all elements and all parameters.We want to impose that or In order to compute a specific function, a fixed point iteration scheme is assumed with all other functions being constants.
In the fixed point scheme only the physical space or one variable are considered at a time, e.g.μ α , so that the test space (the g's) are restricted to the basis used for that variable.
To determine the values of f n for all elements, which define the solution in the physical space, we solve: For a specific function F n α (μ α ) the iteration procedure determines the coefficients F n α that satisfy: Finally, the parametric PGD approximation in terms of the displacements is: and in terms of axial forces: As is typical for the PGD method, the greater the number of terms in the sum, the best should be the approximation achieved.The converge criteria adopted here is the difference between the strain energy of the current pair of solutions (compatible and equilibrium).The process is stopped when either the error is lower than a required tolerance, or when the process stagnates.
When non homogeneous boundary conditions are considered (u(0) = 0 and N (L) = 0, respectively for the compatible and for the equilibrium models) their effect are taken into account in the first term of the sum in (36), so that for the other terms the boundary conditions are homogeneous.

Parametric problem and approximated error measurements
We defined μ as a vector of D parameters μ 1 , μ 2 , . . ., μ D .Since for our specific example the integrals in the physical space are already considered in the system matrices presented in section 2, we can write the parametric finite element approximations of the potential energies directly in terms of the stiffness and flexibility matrices and of the vectors of the applied force and imposed displacement, such that: and We apply the same concept to obtain the bound of the energy of the error, so that: or we can simply define it in terms of the energy density as: When considering the error measurements in the parametric form, we will obviously have a different result for each combination of parameters we use.Therefore, it is necessary to have an additional way to determine the quality of the solutions obtained, that also considers the effects of the parametric domain.One simple solution is to integrate the error measurements obtained in the parametric domain, resulting in a single value that accounts for the quality of all possible solutions.For the energy of the error in particular, it is worth noting that if the integral in space is extended to the parameters, the bounding properties still hold, so that: We can obtain similar approximated error measures using the PGD approximations.We need only to substitute the finite element approximated displacements û and axial forces N for their PGD counterparts ûpgd and Npgd .In the elements domain, a local coordinate system is used (ξ ∈ [−1, 1]), which is linearly mapped to x, which is then linearly mapped to x.

Practical aspects of the discretization of the physical domain
Using φ 0 and φ 1 to represent the linear interpolation functions associated with the end nodes of an interval we have: These transformations are illustrated in Fig. 2. Each section may be divided in several finite elements.Normally, for problems with fixed geometry, only one mapping is used, from the frame of element (ξ ) to the global coordinates (x).
Since both mappings are linear the Jacobian of each transformation is constant, and the Jacobian of the double mapping, J , is equal to their product.The derivative of an arbitrary function, f , is obtained as and its integral is Our basis, in local coordinates, are obtained by combining the linear interpolation functions with the Legendre polynomials L i (ξ ): We can re-write, for example, the stiffness and flexibility matrices in the local frame of the elements, by applying (52) and (53): Fig. 3 Basis for one dimensional approximations The complete system of equations for each problem combines the assembled stiffness (or flexibility) matrices of each element of each section, together with additional constraints that impose displacement continuity (5) or equilibrium of normal forces (9) between sections.This corresponds to impose: • Continuity of displacements and (weak) equilibrium of nodal forces, or • Equilibrium of normal forces and (weak) continuity of the displacements.
For the compatible formulation, to have a continuous displacement at the node between sections b and b + 1, we must impose that: This condition can be expressed in matrix form as The transpose of these matrices reflects the effect of the force at the node (transmitted between sections b and b+1) in the equations of weak equilibrium of the adjacent sections.
For the equilibrated formulation, the process is complementary: we impose continuity of the nodal forces and the displacement of the interface is accounted for in the weak compatibility conditions of each adjacent bar.The matrices involved are the same.
Note that, in case we wanted to consider a concentrated force applied at a node, or a relative imposed displacement, the equations need to be modified accordingly.
In the practical implementation of the problem, we apply these mapping considerations for all the equations described so far, but we will omit them from now on to reduce repetition of the text.

Characterization of the test case
To apply the concepts presented here, we consider the problem shown in Fig. 4, with the following characteristics: • The bar is composed of two sections, with a total summed length L of one; • Each section has its own support stiffness k b ; • The axial stiffness EA for the first section has a unit value; • The axial stiffness of the second section is equal to β;

x1 x2
Fig. 4 A problem with two sections, connected at an intermediate point, defined by γ • The length of the first section is equal to γ .
The parameters for the PGD with their respective limits are: The calculations are performed considering an applied force at the free end (P x=L = 1, x=0 = 0).Unless otherwise stated, the polynomial approximations are linear for all parameters.Notice that these degrees can lead to solutions which are not accurate, but they serve the purpose of this paper, which is to identify regions where the simulations can be improved.
The tolerances τ adopted for the fixed point iteration scheme and the PGD enrichment process are τ fix = 1 × 10 −3 and τ pgd = 1 × 10 −6 , respectively, with the maximum amount of iterations being 3 for the fixed point and 3001 for the PGD.Setting such low number of iterations for the fixed point scheme may lead to more terms for the convergence of the PGD, but drastically decreases the computational time, which is a phenomenon also observed in [11,12], but should not be generalized for all problems, as this behavior is case dependent.
To further improve the convergence of the solutions, the limits for the parameters k 1 , k 2 and β were mapped to a logarithmic scale k 1 = 10 k1 , k 2 = 10 k2 and β = 10 β , with ranges k1 . This mapping does not affect the PGD approximations, serving only to diminish the influence of the edges in the solutions.Additional simulations must be performed with different limits to assess the effect of this modification.
We now present solutions in terms of the error in energy and the energy of the error.The simulations were performed using quadratic approximations in space, in order to better visualize the bounds of the solutions.The bounding characteristics discussed in the previous section are observed in Fig. 5, which presents selected results for our problem.The left part of the Figure shows the integral of the energies of the PGD and FEM approximations, while the right part shows the integral of the energy of the error.The energies for the FEM solutions are very close to the exact solutions and are represented Fig. 5 Integral of the energies for both PGD and FEM approximations as details in the figure.Notice that the exact energy of the error is zero and therefore it is not shown in the figure while the exact integral of the error in energy is greater than zero and is presented in the left figure with a black line.The energies of the compatible and equilibrated solutions will have a convergence from above and bellow, respectively, when the problem is force-driven, which is the behavior observed.On the other hand, for the bound of the energy of the error we have an upper bound, meaning that it should be always greater than the values found for the exact solutions.From the figure it is also possible to see that, as the number of PGD modes increases, the two complementary PGD solutions converge to the approximation of the FEM, always preserving their bounding characteristics.
Figure 6 shows the relative difference between the approximations obtained from the PGD and FEM models.The integral of the energies of the PGD model will converge to the FEM energy as the number of modes in the solution increases.The results obtained for the integral of the energy of the solutions and the energy of the error are virtually the same.This sort of behavior is expected for force or displacement driven problem, due to the orthogonality of the FEM solutions, but is not the case for mixed problems.
Figure 7 shows the behavior of the relative difference between the integral of the energies of the PGD and FEM approximations and the average potential or strain energies as the degree or the number of elements per section used to obtain the solutions increase.The solutions behave as expected, with better results being achieved as the degree or the number of elements increases.

Fig. 6
Relative difference between the PGD and FEM integral of the solutions energies and the average total absolute energy of the FEM solution (left).Relative difference between the PGD and FEM integral of the error energies and the average strain energy of the FEM error (right) 7 Relative difference for the integral of the solutions energy and the integral of the error energy for the PGD approximation with different number of degrees and elements Notice that for higher degree approximations, the solution starts to lose its convergence rate and also the initial solutions are worse than the ones obtained with smaller degrees.The decrease of the convergence rate is caused by a limit in the number of PGD modes due to the tolerance of the simulations.An indication that the solution is being constrained by the tolerance is that the number of PGD modes for a higher degree or number of elements may decrease.This can be seen for the solution with p = 5 and h = 4 and h = 5.Decreasing the tolerance can recover the rate of convergence, although this will not fix the poor results for the initial modes.We believe these poor results are a consequence of the complexity of the higher degree solutions.This is suggested by the fact that this behavior is not observed when increasing the number of elements of the solution.The results in terms of the error in energy and energy of the error are again the same and were omitted from Fig. 7.
The error in energy is simpler to visualize, as it is based in the energy of the solutions which has a physical meaning.Also, each compatible and equilibrated solution has its own energy, allowing for individual analyses of the error.In the other hand, the bounds that are seen in the error in energy are limited to situations when we have force or displacement driven problems, which is not the case for the energy of the error.That being said, from this point on we will base or solutions solely in the energy of the error, as it provides error bounds for general problems.

Solution adaptivity
Adaptivity mesh techniques are designed to improve approximated solutions by modifying either the element mesh configuration or their degrees of freedom [9].One of the key points in the mesh adaptivity process is the proper definition of the relation between the solution convergence and the element size and/or degree of the approximation functions.This convergence dependency is expected to be precise when a smooth solution is being studied and the mesh that represents the solution domain has a large enough number of elements.
We can usually divide the adaptivity process in three categories: the modification of the size of the elements in a mesh (h-adaptivity); the modification of the approximation degree for the functions (p-adaptivity); and a combination of both these process (hp-adaptivity).This paper explores the aspects of the h-and p-adaptivities, but not the combination of both.
We want to obtain an adaptivity indicator that is capable of capturing the best regions to refine, without differentiating between the physical or the parametric space.To achieve this, we collect in variable χ both the parameters and the physical space.This means that, from this point on, a reference to the parameters χ of the problem may be referring to a material parameter μ or to the parameter x of the physical space.
The solutions presented will be shown in terms of the number of elements n h or the degrees of the polynomial approximations p.Our specific problem has a total number of parameters n χ equal to five: four μ's and the physical space, which we divide in the two sections, b 1 and b 2 .

Adaptivity indicator
A question that arises is how do we use the complementary solutions to decide where to refine the model.We know that the integral of the energy of the error covers the whole parameters and physical space domain and therefore it is always the same, no matter in which order we performs the integrals.The idea was, then, to look at the effects of each parameter independently by performing the integral of the error density ρ in all parameters but one.If this function is constant, it means that the parameter is not affecting the error bound, indicating that we need to look at the derivatives of the functions.This lead to the idea of working with the derivative of ρ with respect to all parameters.After numerical tests to assess how to relate the derivative of the augmented set of parameters χ with the effect of that parameter on the error, the following expression for an error indicator ι was selected: where: This expression was selected after testing several dimensionally consistent alternatives, including combinations of the derivatives of ρ with respect to the parameters.From this study, we concluded that (60) is the alternative that best captures the regions where refinement is required.Note that this expression is not used to control if the process has converged, only where to refine next.We have not yet found a theoretical background that supports its application.
Considering ρ, which is non negative, as a pseudo mass density, we can interpret ι(χ) as the first order moment of the derivative of the error with respect to the center of mass of the domain.Recalling that the subscript k indicates solutions that come from kinematically admissible displacements and s that the solutions come from statically admissible axial forces, we can write ρ for an element, considering the domain decomposition particularities, as: where the components in (62) can we approximated in a similar manner as in (28), so that: These approximations can be expressed both in terms of FEM or PGD, only needing to set ûk and Ns accordingly.In order to compute (60) we also need to compute the derivative ∂ρ(ξ ,μ) ∂χ and, therefore, we need the derivatives of (62).We will have different expressions for derivatives depending on what type of parameter χ is representing.
When χ represents the physical domain parameter x, and assuming that the material parameters k and β are defined such that k = μ 1 and β = 1 if x ≤ μ 4 , and k = μ 2 and β = μ 3 if x > μ 4 we have: And when χ represents one of the parameters μ i , with i = 1, 2, . . ., D, we have: Notice that d dμ i dξ dx = 0 unless μ i is the geometric parameter μ 4 = γ .For the derivatives of the PGD approximations, we have: We can now apply Eqs. ( 63) and (65) into Eq.( 60) to obtain the error indicator ι.Notice that by multiplying the derivative by the variable, we are able to keep ι with units that are consistent with the energy of the error 2 .

p-Adaptivity
This section studies the results for p-refinement when using the proposed adaptivity indicator.It works with a simple verification, by identifying the parameter that has the highest value for ι among all the parameters studied and increasing the polynomial degree of that parameter by one, as this parameter is expected to influence the solution error the most.
The adaptivity process using the indicator proposed is compared with the uniform mesh adaptivity and is shown in Fig. 8. Using the adaptivity indicator leads to better results than simply uniformly increasing the degrees of approximation functions, as it chooses the parameter that has the major influence in the solution to improve its degree.The drawback of this process is the need to repeat the simulation for every new degree improved.
One way to determine if the method to drive the adaptivity process is reliable is to determine all the possible solutions for a given set of degrees for the polynomial approximation and verify if the method is capable to accurately define which parameter should have its degree improved.
The following computation process was done: 1 Start with linear functions for the parameters and quadratic for the physical space; 2 Increase the polynomial degree of the parameter by one and compute the error; 3 Return the parameter to the original degree and repeat the previous step for the next parameter, until all the variables are studied; 4 Pick the parameter that leads to the smallest error and permanently increase the polynomial degree of that parameter by one; 5 Go back to step two until the sum of the degrees of the parameter reaches a predetermined value.
The results for different tolerances (τ pgd ) were compared with a reference solution (τ ref ) which was obtained by testing all possible solutions and choosing the optimal result, using a tolerance of 1×10 −12 .The results can be seen in Fig. 9, where the error tends to stagnate after a certain tolerance is achieved.
The adaptivity process is capable of precisely identify the parameter that causes the most impact in the solution, achieving better results for the error sooner than increasing all the degrees at once.We consider that the adaptivity method works if it is capable of selecting the parameter that required its polynomial degree increased to achieve the smallest error.Therefore, as long as the tolerance is small enough, the adaptivity method proposed is capable of choosing the optimal parameter to be refined.

h-Adaptivity
We now study the simulations results for h-refinement when using the proposed adaptivity indicator.Again, just like for the p-refinement, we use the parameter that gives the highest value ι as the indicator mechanism for the adaptivity process.
One key difference for the h-refinement is that it is also necessary to know where the element should be divided.The most natural place to choose is the middle of the element, but we can instead use the center of gravity of the energy error for that parameter, which is already being computed for the adaptivity indicator.To confirm that the center of gravity Fig. 9 Comparison between the reference simulation τ ref and the results with different tolerances for the integral of the energy of the error as a function of the combination of polynomial degrees Fig. 10 Energy of the error for variable β considering different element break positions.Details for the optimal break point and the computed using the gravity center of the energy of the error curve for β of the parameter is the best place to break the element, a simulation with 50 different break points in the interval ] − 1, 1 [ of the variable β was performed and is presented in Fig. 10.It is possible to observe that the point that leads to the real minimum value is really close to the center of gravity of the parameter.This behavior is seen in all others parameters, for different numbers of elements or degrees, leading us to believe that the center of gravity provides a good indication on where to divide the element.
The adaptivity process using the proposed indicator is compared with the uniform mesh adaptivity and is shown in Fig. 11.Using the adaptivity indicator leads to better results than uniformly increasing the number of elements, as it chooses the variable that has the major influence in the solution to divides its elements.Just like for the p-refinement, the drawback of this process is the need to repeat the simulation for every new element division.
We repeat the verification process performed for the p-refinement, determining all possible solutions to the problem for a given sets of elements and verifying if the method is able to accurately define which parameter should have its element divided.The following computation process was performed: 1. Start with all parameter with single element; 2. Divide the element of the parameter and compute the total error; 3.Return the parameter to the original number of elements and repeat the previous step for the next parameter, until all the variables were studied; 4. Pick the parameter that results in the smallest error and permanent split the element of that parameter; 5. Go back to step two until the simulation reaches a predetermined value.
The results for different tolerances (τ pgd ) were compared with a reference solution (τ ref ) which was obtained by testing all possible solutions and choosing the optimal result, using a tolerance of 1 × 10 −9 .These results are presented in Fig. 12, with the total error stagnating after a certain tolerance is achieved.
The adaptivity process is capable of precisely identify the parameter that causes the most impact in the solution, achieving better results for the total error sooner than increasing all the elements at once.We consider the adaptivity method good if it is capable of selecting the parameter that needed its element split to achieve the smallest error.Therefore, just like for the p-refinement, as long as the tolerance is small enough, the adaptivity method proposed is capable of choosing the optimal parameter to be refined.

Conclusion
This paper defines new strategies to assess the error of a PGD parametric problem.We give a brief explanation on how the dual analysis of error works, defining errors measures in terms of the energy of the error and the error in energy.We proceed by introducing finite element and PGD approximations, focused for the solution of a 1D linear elastic problem, considering the details for the implementation of a geometric parameter, specific the length of the sections of the bar we are studying.
We present several examples to evaluate the behavior the PGD approximations.We compute the integrals of the energy of the error and the potential energies in the parametric domain for both PGD and FEM and compare the results with the exact solutions.The integral of the energy of the error presented itself as a good error measure, as it bounds the solutions without any additional conditions, while the error of the energies requires force or displacement driven problems to obtain bounds of the solutions.
We were able to develop an adaptivity indicator, based on empiric data obtained from several simulations performed.This adaptivity indicator is capable of capturing the optimal parameter (or space region) to be refined, both in terms of p-and h-refinements, leading to lower error values than those obtain with a simple uniform refinement.
Additional work is being done to extend the results obtained to a 2D and 3D framework [13].We can also extend the results to obtain quantities of interest and bounds of their errors from for the PGD approximations [14], giving a direct physical meaning to the solutions obtained.

A
global coordinate system (x ∈ [0, L]) is used, where the coordinates of the initial and final nodes of section b are γ b−1 and γ b , with γ 0 = 0 and γ n b = L. Therefore, the geometric parameters correspond to the γ i , with i = 1, . . ., (n b − 1).Each section, which is divided in n e[b] elements, uses an intermediate coordinate system (x [b] ∈ [0, 1]).The coordinates of the initial and final nodes of element e are γ[b]e−1 and γ[b]e , with γ[b]0 = 0 and γ[b]n e = 1.In the following, we replace γ[b]e with γe , x[b]e with xe and x[b] with x, unless the reference to the section of the element is necessary.

Fig. 2
Fig.2Geometric mappings: from bar domain to section and from section domain to element

Fig. 8
Fig. 8 Comparison between the uniform mesh refinement and the adaptivity method for the integral of the energy of the error as a function of the combination of polynomial degrees.The detail numbers are the the degrees of the approximation in each parameter, respectively k 1 , k 2 , β 1 , γ 1 , b 1 , b 2

Fig. 11
Fig. 11 Comparison between the uniform mesh refinement and the adaptivity method for the integral of the energy of the error as a function of the element partitions.The detail numbers are the number of elements in each parameter, respectively k 1 , k 2 , β 1 , γ 1 , b 1 , b 2

Fig. 12
Fig. 12 Comparison between the reference simulation τ ref and the results with different tolerances for the integral of the energy of the error as a function of the element partitions