Denise De

Agnaldo M Farias
email: agnaldo@gmail.com

Phillipe R B Devloo

Mixed finite element approximations of a singular elliptic problem based on

Keywords: Mixed finite elements, hp adaptivity, curved and anisotropic meshes, quarter-point elements, boundary point singularity

some anisotropic and hp-adaptive curved quarter-point elements

Introduction

A source of difficulties in the numerical simulation of elliptic problems in bounded domains is the occurrence of singularities in the solution. For instance, they can be caused by the presence of re-entrant boundary corners (or by abrupt change in boundary conditions), typical situations in fracture mechanics applications, where a radial square root singularity is representative for modeling the state of tension around a fracture in elasticity formulations.

Since the beginning of finite element history, poor convergence of standard H 1conforming methods based on regular grids has been verified in the presence of this kind of low regularity behavior. One popular remedy for this drawback consists in adding special functions to the finite element approximation space, as early described in the papers [START_REF] Fix | On the use of singular functions with the finite element method[END_REF]. This is also the principle guiding the extended finite element method [START_REF] Moës | A finite element method for crack without remeshing[END_REF][START_REF] Daux | Arbitrary branched and intersecting cracks with the extended finite element method[END_REF][START_REF] Kyoungsoo | Integration of singular enrichment functions in the generalized/extended finite element method for threedimensional problems[END_REF][START_REF] Yazid | A state-of-the-art review of the x-fem for computational fracture mechanics[END_REF], using the notion of partition of unity [START_REF] Melenk | The partition of unity finite element method: Basic theory and applications[END_REF]. Another different concept is the application of isoparametric mappings, with singular map for the geometry and for the construction of shape functions, as described in [START_REF] Barsoum | On the use of isoparametric finite element in linear fracture mechanics[END_REF][START_REF] Henshell | Crack tip finite elements are unnecessary[END_REF]. In this context, the quarterpoint element is widely used in linear elastic fracture computations with H 1 -conforming finite element methods. Adaptive numerical treatment of singular problems, incorporating local mesh refinement (h-adaptivity) and/or different orders of approximating polynomials on separate elements in the same mesh (p-adaptivity), has proved to be very effective [START_REF] Ainsworth | A Posteriori Error Estimation in Finite Element Analysis[END_REF][START_REF] Demkowicz | Computing with Hp-Adaptive Finite Elements[END_REF][START_REF] Schwab | hp-FEM for Fluid Flow Simulation[END_REF]. Anisotropic mesh refinement has also been considered and analyzed in [START_REF] Apel | Anisotropic Finite Elements: Local Estimates and Applications[END_REF][START_REF] Farhloul | Some mixed finite element methods on anisotropic meshes[END_REF][START_REF] Acosta | Error estimates for raviart-thomas interpolation of any order on anisotropic tetrahedra[END_REF][START_REF] Li | Regularity and a priori error analysis on anisotropic meshes of a Dirichlet problem in polyhedral domains[END_REF], with improved results. There are also the so called Mapped Finite Element Methods (MFEM), initially proposed in [START_REF] Tsamasphyrosm | The mapped elements for the solution of cracked bodies[END_REF], for better approximations of straight crack problems, and recently extended in [START_REF] Chiaramonte | Mapped finite element methods: High-order approximations of problems on domains with cracks and corners[END_REF] for more general crack configurations. Namely, the method can be interpreted as a reformulation of the problem statement in a reference parametric domain Ω, over which the singularity is eliminated, the corresponding parametric solution û being smoother than the original targeted one, and thus being able to be accurately solved using approximations based on standard uniform partitions of Ω. Similar ideas have been used for the construction of the Method of Auxiliary Mapping (MAM), summarized in [START_REF] Babuška | Finite element method for solving problems with singular solutions[END_REF], which has been applied for p and hp-versions of various elliptic and elasticity problems.

The developed method incorporates all these aspects in some extent, but the adopted formulation is of mixed methods, where two variables are simulated simultaneously: the primal (potential) variable u, approximated by discontinuous functions in L 2 (Ω), and the dual (flux) variable σ, represented in H(div)-conforming spaces [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF]. The interest is to evaluate the performance of this method for the resolution of a square-root radial point singularity by the incorporation of approximations based on quarter-point elements in the vicinity of the singularity point, the shape functions inheriting the singular behavior of the geometric map. As far as we understand, mixed formulations based in quarter-point elements have not been explored in the literature.

Two kinds of quarter-point elements shall be considered: the eight-noded quadrilateral element, placing two mid-side nodes near the singular vertex, and the collapsed quadrilateral quarter-point element, obtained by reducing one of the quadrilateral element edges to the singularity point, resulting a six node triangular element. The square-root radial singularity is exactly captured along element edges by the first type of quarter-point element, while in the latter case, it is captured inside the coarse macro elements as well. Starting from fixed macro quarter-point elements at the coarse level, they are refined by mapping uniformly refined square meshes on the master element by the corresponding geometric transforms. This uniform refinement reveals shape regular curved meshes inside coarse eight-noded quarter-point elements, and anisotropic meshes inside the coarse collapsed quadrilateral coarse elements, grading towards the singular point. In the later case, the meshes are formed by some triangles sharing the point singularity, and trapezoidal quadrilateral elements elsewhere.

As shall be illustrated by some verification tests, the proposed formulation based on quarter-point geometry share some underlying ideas of MFEM and MAM methods, in the sense that the method can be interpreted as local reformulation of the problem statement on the macro elements of the coarse mesh in the reference parametric domain K, over which the singularity is eliminated, the corresponding parametric solution û being smoother than the original targeted one, and thus being able to be accurately solved using approximations based meshes uniformly refined. This effect has a great impact on simulations for the collapsed quadrilateral quarter-point scenario, but it has some limitations for eight-noded quadrilateral ones. For this later case, hp-adaptive meshes are employed for the quarter-point geometry in order to accelerate the convergence of the resulting approximate solutions.

Several methods have been developed for the construction of divergence conforming approximation spaces to be applied in flux approximations of the mixed element formulation. As described in [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Arnold | Quadrilateral H(div) finite elements[END_REF], they are usually formed by locally backtracking vector polynomial spaces restricted to the master element by the Piola transformation, which is defined in terms of the geometric map. The constructions of hierarchical high order spaces in [START_REF] Demkowicz | Polynomial exact sequence and projection-based interpolation with application maxwell equations[END_REF][START_REF] Solin | Higher-Order Finite Element Methods[END_REF][START_REF] Zaglamayr | Hight order finite element methods for electromagnetic field computation[END_REF] are based on the properties of the De Rham complex. For the applications of the present paper, H(div)-conforming finite element subspaces based on curved quadrilateral meshes, with hp-adaptation, are constructed. The vector shape functions form a hierarchy of bases, as proposed in [START_REF] Siqueira | A procedure for the construction of hierarchical high order hdiv and hcurl finite element spaces[END_REF] for triangular and quadrilateral affine elements, and extended to three-dimensional tetrahedral, hexahedral and prismatic affine meshes in [START_REF] Castro | Three dimensional hierarchical mixed finite element approximations with enhanced primal variable accuracy[END_REF]. Cases of hp-adaptive mixed formulations for other kinds of elements in two-and three-dimensional problems, including curved and surface elements, have been treated in [START_REF] Devloo | Two-dimensional hp-adaptive finite element spaces for mixed formulations[END_REF][START_REF] Devloo | Mixed finite element approximations based on 3D hp-adaptive curved meshes with two types of h(div)-conforming spaces[END_REF].

The paper is organized as follows. Some basic facts about the mixed element formulation for elliptic problems are described in Section 2. General aspects on the construction of H(div)-conforming approximation spaces for hp-adapted curved meshes, and the kinds of adopted approximation space configurations, are recalled in Section 3. Some properties of approximations based on quarter-point geometry are illustrated in Section 4. Applications of these space configurations for the mixed formulation of a model problem are discussed in Section 5, for uniform and hp-adaptive frameworks. The results are confronted with the corresponding ones given by simulations based on affine square meshes, showing the dramatic error reduction when quarter-point elements are used. Section 6 gives the final conclusions.

The implementations are performed in the NeoPZ 1 computational platform, which is an open-source object-oriented project providing a comprehensive set of high performance tools for finite element simulations, including the ones required for the present numerical study, namely H(div)-conforming approximation spaces, hp adaptivity, curved elements, and appropriate integration rules for functions showing point singularities.

Mixed finite element formulation

Consider a model Poisson problem:

σ = -∇u in Ω, ∇ • σ = f in Ω, u = u d on Γ D , ∇u • η = g on Γ N ,
where

∂Ω = Γ D ∪ Γ N , Γ D ∪ Γ N = ∅,
u ∈ U = L 2 (Ω) and σ ∈ V = {q ∈ H(div, Ω) : q • η| Γ N = -g}, such that ∀ ϕ ∈ U and ∀ q ∈ H(div, Ω), with q • η| Γ N = 0, a(σ, q) + b(q, u) = c(q), b(σ, ϕ) = (ϕ), (1)
where

a(σ, q) = Ω σ • q dΩ, b(q, u) = - Ω u∇ • q dΩ, c(q) = - Γ d u d q • η ds, (ϕ) = - Ω f ϕ dΩ.
In typical H(div)-conforming discretized versions of this formulation, approximate solutions σ T and u T are searched in finite dimensional spaces V T ⊂ V and U T ⊂ U , and the weak formulation (1) is enforced with test functions in the corresponding spaces ϕ ∈ U T and q ∈ V T , with vanishing normal components over ∂Ω.

The natural functional context for the analysis of the mixed formulation (1), and of its discretized versions, is of saddle-point problems, for which the celebrated inf-sup condition (also known as LBB condition, for Ladyženskaja-Babuška-Brezzi) is necessary and sufficient for well-posedness determination. The verification of the inf-sup condition requires a careful choice of the approximation spaces V T and U T , by enforcing the property ∇ • V = Û of the associated spaces on the master element. This fact allows the definition of a projection π over V commuting the de Rham diagram

H α (Ω) ∇• -→ L 2 (Ω) ↓ π ↓ λ V ∇• -→ Û ( 2
)
where λ is the usual L 2 -projection over Û . These projections are globally extended to Π T over V T , and Λ T over U T .

For such stable finite element approximation spaces, the following estimates hold for the approximate solutions σ T ∈ V T and u T ∈ U T [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Durán | Mixed finite element methods[END_REF]:

||σ -σ T || L 2 ||σ -Π T σ|| L 2 , ||∇ • (σ -σ T )|| L 2 ||∇ • (σ -Π T σ)|| L 2 , ||u -u T || L 2 {||σ -Π T σ|| L 2 + ||u -Λ T u|| L 2 }.
We also remark that, assuming elliptic-regularity (verified by convex polygons or smooth regions Ω), and using duality arguments, enhanced order of accuracy can be obtained for u T when enriched space configurations are used, as in the case of RT + k spaces described next. Projection convergence rates in L 2 -norms for these finite element approximations are usually determined by the largest total degree of polynomials included in the local spaces V(K) = F div K V and Û (K) = F K ( Û ). For instance, assuming shape regular meshes T h , with maximum element diameters h, and uniform polynomial distribution all over the elements, if k and t are the largest natural numbers such that [P k ] 2 ⊂ V(K), and P t ⊂ U (K), then for σ ∈ H q (Ω) and u ∈ H s (Ω), with integer indices 0 < q ≤ k +1, and 0 < s ≤ t + 1, the next estimations hold for approximations based on affine and for bilinearly mapped elements [START_REF] Boffi | Mixed Finite Element Methods and Applications[END_REF][START_REF] Arnold | Quadrilateral H(div) finite elements[END_REF]

||σ -Π h σ|| L 2 (Ω) h q+1 |σ| H k+1 (Ω) , (3)
||u -Λ h || L 2 h t+1 |u| H s+1 (Ω) , (4)
the constants arising in the error estimates being depend on the ratio between the outer and inner diameter of the elements.

In [START_REF] Ern | Finite element quasi-interpolation and best approximation[END_REF], Corollary 5.3, quasi-interpolation operators have been constructed for which similar estimations hold for the cases of fractional Sobolev indices q and s, for affine meshes. Consequently, high order approximation spaces based on affine meshes have limited usefulness for problems with corner sigularities, with slow convergence determined by the low regularity of the solution u ∈ H s (Ω), 0 < s < 1. For such problems, methods using anisotropic mesh grading (with aspect ratio of elements growing to infinity in a refined region close to the singularity) are capable of achieving better accuracy. But for such frameworks, the above mentioned error analysis does not hold. The works in [START_REF] Farhloul | Some mixed finite element methods on anisotropic meshes[END_REF][START_REF] Acosta | Error estimates for raviart-thomas interpolation of any order on anisotropic tetrahedra[END_REF] are dedicated to the study of interpolation properties of Raviart-Thomas elements for some specific anisotropic affine meshes, which are meant for the approximation of this kind of singular solutions. In [START_REF] Farhloul | Some mixed finite element methods on anisotropic meshes[END_REF], error estimates for the mixed formulation for approximations based on anisotropic meshes are obtained, showing first order improvement.

For the current applications using quarter-point elements, the quadratic geometric mappings are not affine nor bilinear. However, the mapped approximation spaces inherit the singular behavior of the geometric map, which is in accordance with the targeted problem solution. As shall be verified, due to the adopted particular refinement procedure, by a global mapping of uniform partitions on the master element, the resulting meshes based on eight-noded quadrilateral quarter-point elements are shape regular, and graded anisotropic meshes occur for collapsed quadrilateral quarter point elements.

Space configurations based on h, p and hp-adaptive curved meshes

Let Ω ⊂ R 2 be a Lipschitz computational domain, with boundary ∂Ω, and denote by η its outward unit normal. Let T be a partition of Ω formed by elements K. We shall be concerned with elements mapped from rectangular master element by quarter-point quadratic geometric transforms, and T does not have, in principle, any limitation on hanging sides and approximation order distribution k = (k K ).

Approximation vector spaces

V k T ⊂ H(div, Ω) = q ∈ L 2 (Ω) 2 ; ∇ • q ∈ L 2 (Ω) ,
which are defined piecewise over the elements of T , require that the local pieces q K = q| K should be assembled by keeping continuous normal components across common element edges. For potential approximations, scalar spaces U k T ⊂ L 2 (Ω) are also piecewise defined over the elements of T , but without any continuity constraint across element interfaces.

Construction of hp-adaptative spaces

The adoted methodology for the assembly of hp-adaptive H(div)-conforming spaces based on curved meshes has been proposed in [START_REF] Devloo | Mixed finite element approximations based on 3D hp-adaptive curved meshes with two types of h(div)-conforming spaces[END_REF] for three dimensional problems. It follows a sequence of steps. 1. To each element K, there is a bijective geometric mapping F K : K → K, associating each point (ξ, η) ∈ K of the rectangular master element K to a point (x, y) = F K (ξ, η) ∈ K. Usually, it is assumed that F K is differentiable, and that the determinant of its Jacobian matrix does not vanish. An isomorphism F K : φ → ϕ, mapping scalar functions of H 1 ( K) to scalar functios of H 1 (K), is induced by the geometric mapping. It also induces a contravariant Piola transformation F div K : q → q, an isomorphism mapping vector fields q ∈ H(div, K) to vector fields q ∈ H(div, K). It is defined by the formula

q = F div K q = F K [ 1 J K DF K q], (5)
where DF K is the Jacobian of the geometric mapping F K , and

J K = det DF K .
There is also an isomorphism

G : μ → µ mapping functions μ ∈ H 1/2 (∂ K) to functions µ ∈ H 1/2 (∂K) such that µ = F K [ 1 Jn μ], J n = J K | t DF -1 K n|, t DF -1 K
being the transpose of the inverse of DF K . Throughout this paper, n denotes the outward unit normal to the faces of an element K ∈ T .

2. Vector polynomial spaces Vk K are defined in K. For the cases to be considered, the parameter k K refers to the maximum degree of the normal components q • n| ∂ K , and to the full polynomial space

[Q k(K),k(K) ] 2 contained in Vk K .
3. A family of hierarchical bases B K k K = { Φ} for Vk K is given, as proposed in [START_REF] Siqueira | A procedure for the construction of hierarchical high order hdiv and hcurl finite element spaces[END_REF]. The principle is to choose appropriate constant vector fields v, based on the geometry of the master element, which are multiplied by an available hierarchic set of H 1 scalar basic functions φ in order to get Φ = φv. There are shape functions of interior type, with vanishing normal components over all element faces. Otherwise, the shape functions are classified of edge type. The normal component of an edge function coincides on the corresponding edge with the associated scalar shape function, and vanishes over the other edges. [START_REF] Kyoungsoo | Integration of singular enrichment functions in the generalized/extended finite element method for threedimensional problems[END_REF]. A hierarchic family of vector bases B K k K = {Φ} is defined over K by the Piola transformation Φ = F div K ( Φ). They inherit the classification of their counterparts on the master element: there are shape functions of interior type, with vanishing normal components over all element edges. Otherwise, Φ is classified as of edge type, its normal component on the associated edge coinciding with the multiplication of the scalar shape function ϕ = F K φ by the normal component of the Piola mapped vector used in its definition, and vanishing over the other edges.

Approximation spaces

V k T ⊂ H(div, Ω) are formed by functions q ∈ [L 2 (Ω)] 2 ,
which are defined piecewise over the elements of T by local functions

q K = q| K ∈ span B K k K ⊂ H(div, K).
They can be assembled to get continuous normal components on the elements interfaces. This property is obtained as a consequence of the particular properties verified by the proposed vector shape functions and of the continuity of the scalar shape functions used in their construction.

Basic facts about the Piola transformation

As already mentioned, the geometric mapping F K induces the contravariant Piola transformation F div K : q → q, mapping vector-valued functions q ∈ H(div, K) to mapping vector-valued functions Φ ∈ H(div, K). Furthermore, if ϕ = F K φ ∈ H 1 (K), the following identities are valid

K q • ∇ϕ dx dy = K q • ∇ φ dξ dη, (6)
K ϕ ∇ • q dx dy = K φ ∇ • q dξ dη, (7)
∂K q • nϕ ds = ∂ K q • n φ dŝ. (8)
Property [START_REF] Barsoum | On the use of isoparametric finite element in linear fracture mechanics[END_REF] says that the divergence of a transformed vector-valued function is given by the the expression

∇ • q = F K 1 J K ∇ • q , (9)
and property [START_REF] Henshell | Crack tip finite elements are unnecessary[END_REF] means that the normal trace Φ • n over ∂K is the image by the mapping G of the normal trace q • n over ∂ K:

q • n| ∂K = 1 J n q • n| ∂ K . ( 10
)
Specific details for quarter-point maps shall be given in Section 4.

Approximation space configurations

Following the developments in [START_REF] Castro | Three dimensional hierarchical mixed finite element approximations with enhanced primal variable accuracy[END_REF][START_REF] Farias | Two dimensional mixed finite element approximations for elliptic problems with enhanced accuracy for the potential and flux divergence[END_REF], we shall consider two stable configurations for approximation spaces to be used for primal u and dual σ variables in discretized versions of the mixed formulation.

• The first case corresponds to the celebrated Raviart-Thomas configuration RT k , where for a given k = k K , the polynomial space on the master element is

VRT k = Q k+1,k × Q k,k+1 . The scalar space ÛRT k = Q k,k is chosen such that the compatibility property ∇ • VRT k = ÛRT k is verified.
• Another type of space configuration, denoted by RT + k , corresponds to an enriched version of RT k . The flux space VRT + k is formed by including in VRT k the internal shape functions of VRT k+1 , while keeping the edge functions with normal components of degree k. The scalar primal approximations are also enriched to

ÛRT + k = Q k+1,k+1 , such that the compatibility property ∇ • VRT + k = ÛRT + k is kept.
For sufficiently smooth solutions, uniform polynomial degree distribution k(k) = k, and shape regular rectangular meshes, it is well know that convergence rates for the mixed finite element approximations based on the RT k space configuration is of order k + 1 for the flux variable, measured by the H(div)-norm, and for L 2 -errors of the potential variable. As shown in [START_REF] Farias | Two dimensional mixed finite element approximations for elliptic problems with enhanced accuracy for the potential and flux divergence[END_REF], for general bilinearly mapped quadrilateral elements, and for enriched RT + k space configurations, the flux error in H(div)-norm is also of order k + 1, and the potential accuracy enhances to order k + 2 when the Ω is convex.

The computational effort in mixed formulations can be reduced by the application of static condensation, the condensed system being solved only for flux edge terms and a constant value for u in each element, all internal fluxes and potential degrees of freedom, excepting a constant value in each element, are condensed. Thus, the dimension of the static condensed matrix coincides for the approximations spaces of types RT k and RT + k .

4 Quarter-point geometry Following [START_REF] Barsoum | On the use of isoparametric finite element in linear fracture mechanics[END_REF], let K = [-1, 1] × [-1, 1] be the master element, and consider the nodes (ξ i , η i ) ∈ K defined in the next Table 1.

i 1 2 3 4 5 6 7 8 coordinates Vertices Mid-side nodes ξ i -1 1 1 -1 0 1 0 -1 η i -1 -1 1 1 -1 0 1 0
Table 1: Vertices and mid-point nodes of the master element K.

Let N i (ξ, η) be the quadratic Lagrange shape functions

N i (ξ, η) = (1 + ξξ i )(1 + ηη i ) -(1 -ξ 2 )(1 + ηη i ) -(1 -η 2 )(1 + ξξ i ) ξ 2 i η 2 i 4 +(1 -ξ 2 )(1 + ηη i )(1 -ξ 2 i ) η 2 i 2 + (1 -η 2 )(1 + ξξ i )(1 -η 2 i ) ξ 2 i 2 ,
at the nodes (ξ i , η i ). Precisely, for the vertices (ξ i , η i ), where

ξ 2 i = η 2 i = 1, N i (ξ, η) =        1 4 (1 -ξ)(1 -η) [-1 -ξ -η] , at node 1, 1 4 (1 + ξ)(1 -η) [-1 + ξ -η] , at node 2, 1 4 (1 + ξ)(1 + η) [-1 + ξ + η] , at node 3, 1 4 (1 -ξ)(1 + η) [-1 -ξ + η] , at node 4. (11)
For the nodes 5 and 7, with ξ 2 i = 0 and η 2 i = 1,

N i (ξ, η) = 1 2 (1 -ξ 2 )(1 -η), at node 5, 1 2 (1 -ξ 2 )(1 + η), at node 7. ( 12
)
For the nodes 6 and 8, with ξ 2 i = 1 and η 2 i = 0,

N i (ξ, η) = 1 2 (1 -η 2 )(1 + ξ), at node 6, 1 2 (1 -η 2 )(1 -ξ), at node 8. ( 13
)
Given a quadrilateral element K, with vertices (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ), and (x 4 , y 4 ), consider mid-side nodes (x 5 , y 5 ), (x 6 , y 6 ), (x 7 , y 7 ), and (x 8 , y 8 ) placed on the edges (x 1 , y 1 )-(x 2 , y 2 ), (x 2 , y 2 )-(x 3 , y 3 ), (x 3 , y 3 )-(x 4 , y 4 ), and (x 4 , y 4 )-(x 1 , y 1 ), respectively. Define F K : K → K by the formulae

x(ξ, η) = 8 i=1 N i (ξ, η)x i , y(ξ, η) = 8 i=1 N i (ξ, η)y i , such that (x i , y i ) = F K (ξ i , η i ).
According to [START_REF] Barsoum | On the use of isoparametric finite element in linear fracture mechanics[END_REF], the name quarter-point element comes when two mid-point nodes, namely, (x 5 , y 5 ) and (x 8 , y 8 ), of the sides sharing a vertex, (x 1 , y 1 ) for this case, are moved towards it, at one fourth of their associated side length, while the other two mid-side nodes ((x 6 , y 6 ) and (x 7 , y 7 )) are the middle points of the corresponding sides.

Uniform subdivision of a quarter-point element

Given a quarter-point element K mapped by the quadratic transformation F K , its refinement into four sub-elements K j results after the following steps:

1. Let S j be the four squares obtained by a uniform subdivision of K, and consider affine transformations

L j : K → Ŝj 2. The sub-elements in K correspond to K j = F K Ŝj .
3. The quadratic mappings F K j : K → K j associated to K j are defined by the composition

F K j = F K • L j .
Therefore, the vertices of the resulting 4 curved sub-elements in K are determined by the four vertices of K, namely (x 1 , y 1 ), (x 2 , y 2 ), (x 3 , y 3 ), and (x 4 , y 4 ), the four side nodes (x 5 , y 5 ), (x 6 , y 6 ), (x 7 , y 7 ), and (x 8 , y 8 ), and the central point (x 9 , y 9 ) = F K (0, 0). Note that

N i (0, 0) = - ξ 2 i η 2 i 4 + (1 -ξ 2 i ) η 2 i 2 + (1 -η 2 i ) ξ 2 i 2 ,
giving N i (0, 0) = -0.25 for the vertices i = 1, 2, 3 and 4, and N i (0, 0) = 0.5, for the mid-side nodes i = 5, 6, 7 and 8. Consequently

x 9 = -0.25 4 i=1 x i + 0.5 8 i=5 x i , y 9 = -0.25 4 i=1 y i + 0.5 8 i=5 y i .

Example of quarter-point element with 8 nodes

Suppose that K = [0, 0.5] × [0, 0.5], with vertices and mid-side nodes given in Table 2, and illustrated in Figure 1. For this element, the quarter-point map becomes

x = 1 16 (3 + 4ξ + η + ξ 2 -ξ 2 η), y = 1 16 (3 + ξ + 4η + η 2 -ξη 2 ).
Observe that, along the line η = -1, the result of the the quarter-point map gives x = 1 8 (1 + ξ) 2 , and y = 0. That is, over the side y = 0, the square root singular function √ 8x = 1 + ξ is reproduced in the master element side η = -1 by a first order polynomial. Take the more general square-root function w(r, θ) = √ r = (x 2 + y 2 )

1 4
restricted to K, as plotted in Figure 2 (left side). Its representation ŵ(ξ, η) back in the master element K, i.e., ŵ(ξ, η) = w(F K (ξ, η)), has a Taylor expansion around (-1, -1) is given by

ŵ(ξ, η) = η + 1 2 √ 2 + 1 4 √ 2 - η + 1 8 √ 2 (ξ + 1) + 1 16 √ 2 (ξ + 1) 2 (η + 1) + O(η + 1) 3 + O(ξ + 1) 3 + O((η + 1) 3 (ξ + 1)) + O((η + 1) 3 (ξ + 1) 2 ).
Its plot is shown in Figure 2 (right side), which looks much smoother than w, illustrating the capability of quarter-point map F K to well represent square-root type of singularities.

Figure 2: Example: w(r, θ) = √ r, restricted to K (left side), and its parametric representation ŵ(ξ, η) in the master element (ξ, η) ∈ K.

Piola transformation

The Jacobian matrix of F K for this particular quarter-point element is

DF K = 1 4 + 1 8 ξ(1 -η) 1 16 + 1 16 η 2 1 16 (1 -η 2 ) 1 4 + 1 8 η(1 -ξ) .
Note that, over the side y = 0, ∂x ∂ξ = 1 4 (1 + ξ) = 1 4 √ 8x, making the Jacobian singular at the vertex (x 1 , y 1 ) = (0, 0).

By the Piola transformation q = F div K q, the parametrized version q associated to a vector function q has the formula

q(ξ, η) = J K (x)DF K (ξ, η) -1 q = ∂y ∂η -∂x ∂η -∂y ∂ξ ∂x ∂ξ q.
Some examples of parametrized vector functions are shown in Table 3 for this particular eight-noded quadrilateral quarter point element K. Observe that, for q = [1, 0] t and q

= [0, 1] t , q ∈ VRT 1 = Q(2, 1) × Q(1, 2), implying that V RT 1 (K) contains the constant vectors. q q Quadrilateral Q-P Collapsed quadrilateral Q-P 1 0 1 4 + 1 8 η(1 -ξ) -1 16 (1 -η 2 ) 1 8 (1 + ξ) 2 -1 4 η(1 + ξ) 0 1 -1 16 + 1 16 ξ 2 1 4 + 1 8 ξ(1 -η) 0 1 4 (ξ + 1)
Table 3: Parametrized vector functions q associated to constant vector functions q defined in quadrilateral eight-noded and collapsed quadrilateral six-noded quarter-point elements K.

Uniform refinement of K

Setting T K 0 = K, a recursive application of the refinement procedure generates curved refined meshes T K j of the K obtained as the image of the uniformly refined mesh T K j , with spacing 2 -j , by the global quarter-point mapping F K of the macro-element K. The partition T K 1 of the first subdivision step for this particular example of quarterpoint element K is illustrated on the right side of Figure 1, whose middle point is (x 9 , y 9 ) = F K (0, 0) = (0.1875, 0.1875). Illustrations of T K j , for j = 2, 3, are shown in Figure 3 (left side). It can be observed that the uniform refinement of the quarter-point element K produces a directional refinement pattern pointing towards the singularity point (x 1 , x 1 ). Despite of their stretched aspect, the meshes T K j are shape-regular, as illustrated in Figure 6 (left side).

Collapsed quadrilateral quarter-point element

Consider the node numbering on the master element in Table 4, as in [START_REF] Barsoum | Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements[END_REF].

i 1 3 5 7 2 4 6 8 coordinates Vertices Mid-points ξ i -1 1 1 -1 0 1 0 -1 η i -1 -1 1 1 -1 0 1 0 Table 4:
Numbering of vertices and mid-point nodes of the master element K [START_REF] Barsoum | Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements[END_REF].

Suppose that K is the triangular element with six nodes shown in Table 5, which can be interpreted by the collapse of the nodes 1 and 7 in the master element on the singular node 8. For this element, the quarter-point map (x, y) = F K (ξ, η), and the Jacobian matrix of

F K become x = 1 8 (1 + ξ) 2 , y = 1 8 η(1 + ξ) 2 , DF K =    ∂x ∂ξ ∂y ∂ξ ∂x ∂η ∂y ∂η    =   1 4 (1 + ξ) 1 2 η (1 + ξ) 0 1 8 (1 + ξ) 2   ,
with determinant J K = 1 32 (1+ξ) 3 . Parametrized vector functions q associated to some examples of vector functions q defined in K are displayed in Table 3, also showing that V RT 1 (K) contains the constant vectors.

The parametric representation of the radial square-root function w = r 1/2 restricted to K is the infinitely smooth function

ŵ(ξ, η) = 1 2 √ 2 (1 + ξ)(1 + η 2 ) 1/4 ,
which is plotted in Figure 4 (left side). For q = ∇w, restricted to K, the corresponding parametric representation q is given by the regular vector

q = 1 4 √ 2 (1+ξ) (1+η 2 ) 3 4
0 .

The first component of q is represented in the Figure 4 (right side). These facts make the collapsed quadrilateral quarter-point elements more efficient for the representation of radial square-root singularities of type √ r than the quadrilateral eight-noded quarter-point elements, as observed in [START_REF] Barsoum | Triangular quarter-point elements as elastic and perfectly-plastic crack tip elements[END_REF] for linear elasticity, and illustrated by the numerical simulations of the next section. The first and twice subdivisions of the collapsed quadrilateral quarter-point element K are shown in Figure 5. Since x does not depend on η, and fixing constant the value of ξ, the mapped curve is a vertical line. Otherwise, for constant value of η, the mapped curve is the strait line y = ηx. Consequently, the resulting meshes T K j , obtained by applying j times the uniform refinement procedure to the original collapsed quadrilateral quarter-point element K, are formed by j triangles touching the singular point (0, 0), and trapezoidal quadrilateral elements elsewhere. They are graded anisotropic meshes, with growing aspect ratio, as shown in Figure 6 (right side). Figure 6: Aspect ration of the quarter-point meshes T K j obtained by uniform quadrilateral element K, with 8 nodes (left side), and the collapsed quadrilateral quarter point element K, with 6 nodes.

Numerical simulations

A model problem is considered with f = 0, and Ω = [-0.5, 0.5] × [0, 0.5], where the exact solution, shown in Figure 7, is the harmonic function

u = 2 1/4 √ r cos( θ 2 ) = 2 -1/4 x + x 2 + y 2 ,
that presents a radial square root singularity at the boundary point (x, y) = (0, 0) (r = 0), where there is a change from Dirichlet boundary condition u = 0, for x < 0, to Neumann condition ∇u • η = 0, for x > 0. Elsewhere, Neumann boundary condition is enforced. This problem has been used in [START_REF] Forti | A comparative numerical study of different finite element formulations for 2D model elliptic problems: continuous and discontinuous galerkin, mixed and hybrid methods[END_REF] to compare the performance of different finite element formulations, including continuous, discontinuous, mixed and primal hybrid finite element methods. As expected, the application of hp-refinement improves considerably the performance of all methods. The purpose here is to evaluate the efficiency of the mixed method applied to this singular problem when quarter-point elements are used in the vicinity of the singularity.

Numerical results with uniform refinement

Two scenarios shall be considered for the generation of sequences of meshes T j , j = 0, 1, 2, • • • , 5, by using uniform refinement procedures, as described in Section 4.1. Starting from fixed macro quarter-point elements at the coarse level, they are refined by mapping uniformly refined square meshes on the master element by the corresponding geometric transforms.

Mesh scenario 1

The coarsest mesh T 0 is composed of two square elements K 0 , sharing the singular point (0, 0) as one of their vertices, and eight-noded quarter-point maps are applied, as described in the example in Section 4.2. The illustration of the mesh T 2 , obtained by refining twice these coarse elements, is shown in Figure 8 (left side). As already described in the previous section, the resulting sequence of meshes T j is obtained by mapping uniformly refined square meshes on the master element by the corresponding geometric transform of the corresponding father elements K 0 . Consequently, T j is shape-regular. Furthermore, by their definition, the shape functions are produced by the quarter-point geometric map are present in all the elements. However, only two elements neighboring the singular node (0, 0) have a singular mapping. Having in mind the example discussed in Section 4.2, the exact singular solution u = u| K 0 can be represented back in the master element by smoother functions û (ξ, η) = u (F K 0 (ξ, η)) obtained by composing u with the quarter-point mappings of the two coarse elements F K 0 . Consequently, the application of the mixed method to solve the model test problem using stable approximation spaces based on T j are expected to give better convergence rates than using full square uniform meshes. As remarked in the introduction, this behavior is in accordance with the underlying idea of mapped element methods [START_REF] Tsamasphyrosm | The mapped elements for the solution of cracked bodies[END_REF][START_REF] Chiaramonte | Mapped finite element methods: High-order approximations of problems on domains with cracks and corners[END_REF], in the context of H 1 conforming methods for linear elasticity of crack problems, by interpreting the method as a reparametrization of the problem that removes the singularity, the corresponding smoother parametric solution û, with lower gradients, being able to be accurately solved using approximations based on standard affine uniform partitions in the parametric domain.

Figure 9 presents L 2 -error curves for the approximations of u and σ obtained with the mixed method based on RT k and RT + k space configurations based on the meshes T j , 1 ≤ j ≤ 5, with uniform polynomial degree distribution k = 1, 2 (indicated by QP RT k or QP RT + k ). The mesh size h j are for the associated uniform meshes on the master element, whose images by F K 0 form the refined meshes T j restricted to each macro element K 0 . Results for uniform square meshes are also included (denoted by Square RT k ). As expected for this kind of singular solution u ∈ H of flux and potential approximations given by the application of the RT k space configuration based on uniform square meshes, are limited to rates of order 0.5 and 1.0, respectively, independently of the polynomial degree k ≥ 1. When quarter-point elements are applied, enhanced approximations are observed, not only by increasing the rates to order 1 for the flux and 2 for the potential, but also by the significant reduction in error magnitudes. The enriched RT + k space configurations based on quarter-point meshes give more accurate solutions, specially for the potential variable. For instance, the errors in u using RT + 2 on the most refined mesh T 5 is about four orders of magnitude less than using RT 2 on the square mesh with the same resolution.

Figure 8 also includes plots for the results obtained with RT 2 space configuration based on meshes where only the two refined quadrilateral elements sharing the singular vertex are mapped by the quarter-point transformation (indicated by QP * RT 2 ). Confronted with the simulation for RT 2 spaces on the full affine square meshes, the introduction of these two quarter-point elements significantly reduces the error magnitudes, but the convergence rates do not improve. This behavior is in accordance with the directions given in [START_REF] Chiaramonte | Mapped finite element methods: High-order approximations of problems on domains with cracks and corners[END_REF], where the singularities are localized in a fixed sector, where they are properly parametrized, and approximated by uniform discretizations.

Mesh scenario 2

A second mesh scenario considers the coarsest mesh T 0 composed of 4 triangular elements K 0 sharing the singularity point (0, 0), and six-noded quarter-point maps are applied to them. The meshes T j are constructed by the uniform refinement procedure applied to K 0 , as described in the example in Section 4.3. The resulting anisotropic mesh T 2 , obtained by refining twice the coarsest mesh, is shown in Figure 8 (right side).

As in the previous scenario, the exact singular solution u = u| K 0 , restricted to each of the four coarse triangular elements, can be represented back in the master element by smoother functions û (ξ, η) = u (F K 0 (ξ, η)) obtained by composing u with the corresponding quarter-point mappings. However, opposed to the mesh scenario 1, having a limited smoothness improvement, û (ξ, η) are infinitely smooth in the mesh scenario 2, as has been illustrated in Section 4.3. Consequently, the application of the mixed method to solve the model test problem using stable approximation spaces based on the meshes T j of scenario 2 are expected to give better convergence rates than using the meshes of scenario 1. Indeed, as confirmed by the error plots of Figure 10, this property is crucial to compensate the effect of mesh anisotropy, and to produce convergence rates of order k +1 for flux σ and potential u for RT k space configurations, which is typical of standard regular scenarios (shape-regular affine quadrilateral meshes and smooth solutions), and enhanced potential order of convergence k + 2 occurring for the enriched space configuration RT + k . In order to confirm this tendency to higher polynomial degrees, the results in Figure 11 are for a p-refinement experiment based on the coarse level T 1 , clearly showing an exponential error decay as k increases.

hp-Adaptive simulations

Even though the singular solution being better represented by approximations based on quarter-point elements, for uniformly refined partitions T j , as discussed in Section 5.1, Flux error the experiments for the mesh scenario 1, using RT k and RT + k space configurations using eight-noded quadrilateral geometric maps, the improvements revealed to be limited to order k = 1 for the flux, and k = 2 for the potential variables, already verified for k = 1. Now the purpose is to explore an hp-adaptive strategy, combined with these kind of quarter-point elements to enhance the convergence of the approximate solutions. The construction of the hp-adaptive meshes is illustrated in Figure 12. Initially, there is a mesh composed of 8 quadrilateral elements with uniform degree k = 2 distribution. The two quadrilaterals near the singularity are divided into 4 sub-elements and their approximation order are set to k = 3. Then, the refinement process follows the steps:

1. All elements have their approximation order increased by 1.

2. The elements near the singularity are divided and their approximation order are increased by 1.

Two cases shall be compared: one with affine mapping applied to all elements, and the other one for eight-noded quarter-point elements in the central region Ω C . Our purpose is to use this kind of adaptive meshes, for the simulation of the test problem by the mixed formulation using RT k and RT + k space configurations, where k = (k K ) refers to the non-uniform polynomial degree distribution used in the construction of the approximation spaces.

Figure 13 shows the calculated L 2 -norms of the dual and primal errors using these sequences of hp-adaptive quarter-point curved meshes (continuous lines) versus the number of equations solved after static condensation, using RT k and RT + k (o) space configurations. For comparison, results for similar square hp-meshes (dashed lines), without using quarter-point elements, and for uniform meshes with constant k = 2 distribution (+) are plotted. The results demonstrate exponential rates of convergence when using hp-adaptive meshes. However, for the standard hp-adaptive method, based on square elements, the results for RT k and RT + k (o) space configurations are quite the same. Furthermore, they only outperforms the uniform setting with quarter-point elements and constant polynomial degree k = 2 at the very refined levels. Dramatic error reduction is observed when the hp-adaptive strategy is combined with quarterpoint elements in the central region, the accuracy in the primal variable being further improved when the RT + k configuration is applied. Results for square hp-meshes, but inverting the polynomial degree distribution ( ) are also plotted. This type of degree grading, using low order k = 2 on the two elements sharing the singularity point, and increasing the degree as the elements are placed away from the singularity, has been suggested in [START_REF] Schwab | hp-FEM for Fluid Flow Simulation[END_REF], Section 5.2, but for the present model problem the results are less accurate than using the proposed opposite degree grading.

The effect of static condensation is also verified for size reduction of the global system to be solved, being more significant with increasing order of approximation. At the finest level of mesh refinement, the number of condensed equations amounts to more than 94% of the total number of equations.

Final conclusions

This article presents an application of H(div) mixed approximations, combined with singular geometric map applied to the approximation a problem exhibiting singular behavior caused by a abrupt change in boundary conditions.

It was shown that the quality of H(div) approximations of singular problems can be improved dramatically by the use of well chosen geometric maps. This result is an extension to H(div) conforming spaces of known results when using H 1 approximation spaces.

When approximating a problem with square root singularity, the well-known quarter point geometric map leads to an increase in the convergence rate of pressure and flux variables, from order 1 and 0.5 when using uniform meshes, to order 2, and 1, respectively. It is shown that the combination of quarter point geometric maps with h -p adaptive meshes leads to exponential convergence rates.

Even better approximation properties are obtained by using collapsed quadrilateral elements combined with quadratic geometric maps with quarter point singularity. In this case it is demonstrated that the solution in the inverse image of the geometric map is regular. This property is reflected on the accuracy of the numerical results for flux and pressure: rates of h-convergence equivalent to standard regular contexts, of smooth solutions discretized on uniform affine meshes, and exponential rates of convergence when applying uniform p refinements.

It is probable that the results obtained can be extended to three dimensional simulations of problems with singularities.

Figure 1 :

 1 Figure1: Example of quarter-poin-element K: mid-side node movement towards the vertex (x 1 , y 1 ) to one quarter of the correspondind side length l = 0.5 (left side), and the result of one subdivision of K (rigth side).

Figure 3 :

 3 Figure 3: Meshes T K j , j = 2, 3, obtained by uniform subdivisions of eight-noded quarter point element K.

Figure 4 :

 4 Figure 4: Parametric representation ŵ(ξ, η), -1 ≤ ξ ≤ 1, 0 ≤ η ≤ 1, of w(r, θ) = r 1/2restricted to the collapsed quadrilateral six noded quarter point element K (left side), and of the first component of q associated to q = ∇u (the second gradient component vanishess) (right side).

Figure 5 :

 5 Figure 5: Meshes T K j , j = 1, 2, obtained by one and two uniform subdivisions of the collapsed quadrilateral quarter-point element K.

Figure 7 :

 7 Figure 7: Exact solution u = 2 1/4 √ r cos( θ 2 ).

Figure 8 :

 8 Figure 8: Mesh scenarios: T 2 starting with two eight-noded quarter-point elements sharing the singularity point (left side), or starting with four triangular six-noded quarter-point elements.

3 2 - 2 +Figure 9 :

 229 Figure 9: L 2 -erros in u (left side) and in σ (right side) using Raviart-Thomas spaces RT k , or RT + k , with k = 1, 2, based on square meshes (S-dashed lines), on meshes formed by quadrilateral quarter-point elements of 8 nodes (QP-continuous lines), or with just two quadrilateral quarter-point elements of 8 nodes sharing the singular vertex (QP*).

Figure 11 :

 11 Figure 11: L 2 -errors in u (dashed line (s)) and in σ (continuous line (o)) versus the number of equations of the condensed systems, using Raviart-Thomas spaces RT k , k = 1, • • • , 9, based on the fixed collapsed quadrilateral quarter-point mesh T 1 .

Figure 12 :

 12 Figure 12: Illustration of the hp meshes. Affine meshes (left side), and quarter-point meshes (right side), at refinement levels 1, 2, 5, and 6 (from top to bottom).

Figure 13 :

 13 Figure13: L 2 -errors in u (left side), and in σ(right side), using curved quadrilateral meshes (QP) with quarter-point elements (continuous lines), and square (S) meshes (dashed lines), for the mixed formulation based on hp-adaptive RT k , and RT + k (circle) space configurations. Uniform simulations (+) are for square and rectangular quarterpoint elements with RT 2 space distribution. The hp-simulations with inverse degree grading ( ) are for square meshes (SI).

  Γ D and Γ N denoting the parts of ∂Ω for the enforcement of Dirichlet and Neumann boundary conditions, respectively, and η is the outward unit normal to ∂Ω.This problem can be expressed in the form: to find

Table 2 :

 2 Example: vertices and mid-point nodes of a eight-noded quadrilateral quarterpoint element K.

	1/4 l	3/4 l

Table 5 :

 5 Example: vertices and mid-point nodes of a six-noded collapsed quadrilateral quarter-point element K.

  Figure 10: L 2 -errors in u (left side) and in σ (right side), versus the mesh width h of the master element, using Raviart-Thomas spaces RT k , or RT + k , based on quarter-point elements: refining the two quarter-point elements with 8 nodes (QP Q ) (continuous lines), or the four colappsed quadrilateral quarter-point elements with 6 nodes (QP T r ).

				1
			1
	QP Q	RT 1
	QP Q QP Q	+ RT 1 RT 2	2	1
	QP Q	RT 2 +
	QP Tr RT 1
	QP Tr	RT 1 +		1
	QP Tr	RT 2	3
	QP Tr	RT 2 +
	QP Tr	RT 3
	QP Tr	RT 3 +		1
			4

http://github.com/labmec/neopz

Acknowledgements

The authors P R B Devloo and S M Gomes thankfully acknowledge financial support from FAPESP -Fundação de Apoio à Pesquisa do Estado de São Paulo, Brazil (grant 2016/05144-0), and from CNPq -Conselho Nacional de Desenvolvimento Cientítico e Tecnológico, Brazil (grants 305425/2013-7, 305823-2017-5, and 304029/2013-0, 306167/2017-4). P R B. Devloo also acknowledges financial support from ANP-Brazilian National Agency of Petroleum, Natural Gas and Biofuels (grant 2014/00090-2).


Martha Liliana Torres-Barreto

Marianela Luzardo Briceño
email: marianela.luzardo@upb.edu.co

Mileidy Alvarez-Melgarejo
email: mileidyalvarez@hotmail.esa003

Keywords: Capabilities, innovation, public funding, resources Capabilities, innovation, public funding, resources. Resumen

The R&D funding scenario: Can capabilities foster product innovation of firms?

de muchas agendas como Strategic European Research (European Commision, 2013). En este artículo se desarrolla un modelo analítico del impacto tecnológico del financiamiento público para actividades privadas de I + D en términos de innovaciones de producto logradas, para estudiar si sus efectos diferenciales entre empresas dependen en cierta medida de los recursos y capacidades de la empresa. Este análisis se basa en los datos a nivel de empresa de la industria de la manufactura bajo ESSE (encuesta española para estrategia de las empresas), para el período 2000-2012. Se identificaron las actividades tecnológicas de 3350 empresas españolas, se caracterizaron sus recursos y capacidades e identificaron las más relevantes en cuanto a su influencia en la innovación de productos. Se encontró que tanto los fondos públicos de I + D como la inversión privada ayudan en la generación de nuevos productos; pero también, que las capacidades de innovación y la cooperación interempresarial afectan positivamente el resultado tecnológico.

Palabras clave: Capacidades; innovación; financiamiento público; recursos.

INTRODUCTION

The financing of R&D activities is being increasingly examined as public bodies continuously demand measures of impact obtained from the public funding of private R&D activities [START_REF] Giebe | How to allocate R&D (and other) subsidies: An experimentally tested policy recommendation[END_REF][START_REF] Heijs | Freerider behaviour and the public finance of R&D activities in enterprises: The case of the Spanish low interest credits for R&D[END_REF]. To quantify such an impact, a set of economic, academic and technological measures have been proposed. From the technological perspective, the product innovation it's been one of these measures.

Product innovations encompasses the implementation phase [START_REF] Damanpour | Organizational Innovation: a Meta-Analysis of Effects of Determinants and Moderators[END_REF][START_REF] Wolfe | Is firm size conducive to R&D choice? A strategic analysis of product and process innovations[END_REF] and refers to the development and launching of products which are unique or distinctive from existing ones [START_REF] Higgins | Innovate or Evaporate: Creative Techniques for Strategists[END_REF]OECD, 2005). In competitive conditions, they are often considered a crucial resource since the business environment is characterized by rapid technological change, shortened product life cycles and globalization, forcing firms to be more innovative in their product introduction if they want to gain competitive advantages [START_REF] Argüello | Uso de las TIC en las PYME del sector industrial del Área Metropolitana de Bucaramanga[END_REF][START_REF] Chen | IT capabilities and product innovation performance: The roles of corporate entrepreneurship and competitive intensity[END_REF][START_REF] Muralla | Revisitando estrategias de sostenibilidad de las empresas a través de una visión sistémica empresarial[END_REF][START_REF] Song | The Determinants of Japanese New Product Successes[END_REF][START_REF] Torres-Barreto | Exploring the boosting potential of intellectual resources and capabilities on firm´s competitiveness[END_REF]. Their importance is consequently based on the possibilities of entering new markets and take advantages of new opportunities in order to make greater profits.

The link between R&D investments and product innovations has been explored in the literature so far. R&D projects have been identified as precursors of new products [START_REF] Balachandra | Factors for success in R&D projects and new product innovation: a contextual framework[END_REF], while public subsidies are considered to influence decisions of companies in order to improve their innovation processes [START_REF] Afcha | Public funding of R&D and its effect on the composition of business R&D expenditure[END_REF]. Public funding fostering private R&D activities has been identified to positively affect firm innovations through the complementarity effect [START_REF] Griliches | R&D and Productivity: Measurement Issues and Econometric Results[END_REF]; until the point that it can be perceived as a signal of quality for a certain project [START_REF] Beaudry | Impact of public and private research funding on scientific production: The case of nanotechnology[END_REF]. However, succeed in innovating is a complex and difficult task, and it is not easy to explain why few R&D projects come up with innovations that reach the market, while the rest of them fail. Since the task seems to be of a certain difficulty level [START_REF] Eisenhardt | Dynamic capabilities: what are they?[END_REF], we draw on the literature to examine the causes that may lead to some firms to succeed in product innovations and found that capabilities of firms [START_REF] Barney | Firm Resources and Sustained Competitive Advantage[END_REF], in some extent, may contribute to foster the generation of product innovations. Previous results show a positive influence of: cooperation among partners (Antolín-López, Céspedes-Lorente, García- [START_REF] Antolín-López | Fostering product innovation: Differences between new ventures and established firms[END_REF][START_REF] Un | R&D collaborations and product innovation[END_REF], the capability to implement management control systems [START_REF] Bisbe | The effects of the interactive use of management control systems on product innovation[END_REF], capabilities of firms to generate organizational innovations [START_REF] Ballot | The fateful triangle: Complementarities in performance between product, process and organizational innovation in France and the UK[END_REF] and firm´s IT capabilities [START_REF] Bharadwaj | A Resource-Based Perspective on Information Technology Capability and Firm Performance: An Empirical Investigation[END_REF][START_REF] Chen | IT capabilities and product innovation performance: The roles of corporate entrepreneurship and competitive intensity[END_REF]. Consequently, this research attempts to identify the relationship between public funding for R&D activities and the obtention of product innovations, as well as identify capabilities leading to successful product innovations when funding is given to firms.

Using data from the Spanish Survey of Business Strategy (years: 1998Strategy (years: -2013)), we explore the effects on product innovation in the presence of different combinations of capabilities: innovation planning, updating the R&D knowledge base, R&D cooperation and absorptive capability.

With this research, we aim to contribute to theory and practice in several ways. Since empirical studies of the effects of public funding policies on firm competitiveness and innovative output remains limited [START_REF] Chen | IT capabilities and product innovation performance: The roles of corporate entrepreneurship and competitive intensity[END_REF]Afcha & Leon Lopez 2014), and the validity of patents when measuring public R&D funding is still a gap in the literature (Zvi [START_REF] Griliches | Patent Statistics as Economic Indicators: A Survey[END_REF][START_REF] Torres-Barreto | Technological impact of R&D grants on utility models[END_REF], we study certain firm characteristics that could make R&D funding profitable, discuss them as possible drivers in the obtention of innovations that reach the market, and provide a framework for future research. By examining these relationships, we support manager's performance since they need to know whether R&D funding may function as an antecedent for product innovation and if so, under which conditions, so that they can ensure those conditions exist when the funding is received. Moreover, policy makers would be able to design programs that maximize R&D public budget for private research. Is therefore important to explore how this variable interact and what factors may influence this relationship.

THEORETICAL FRAMEWORK AND HYPOTHESIS

PRODUCT INNOVATIONS AND R&D EFFORT

The importance of studying the impact of public R&D funding relies not only on the amount of public money spent but also on the intended effects in both: firms and the society, as product innovation can enhance firms' competitive advantage; the offspring of products changes the technological ecosystem, and, new product introduction redefines patterns, habits and behaviors in society in general [START_REF] Ahlstrom | Innovation and Growth: How Business Contributes to Society[END_REF].

Moreover, product innovations as a technological indicator of R&D activity are also appropriate if the objective is to measure the effect of SMEs´s R&D effort. The process innovation behavior of big incumbent firms contrasts with that of small and medium enterprises who attempt to differentiate via new products [START_REF] Malerba | Schumpeterian patterns of innovation are technologyspecific[END_REF]. Due to this heterogeneity in the innovative activities of firms, big companies are still dominant on post-innovation markets, while smaller ones are likely to be pioneers on the new-product markets (Yin & Zuscovitch, 1998). Furthermore, product innovation can eventually stimulate private R&D investment, which would contribute to minimize the market failure common in the appropriation of R&D [START_REF] Callejón | Public subsidies to business R&D: Do they stimulate private expenditures?[END_REF].

In the quest for measuring the technological impact of public R&D investment, researchers can rely on outcomes as the function of new or improved skills, by a technological gap closed, by the use of a newly developed technology, or by the subsequent technological progress achieved [START_REF] Sohn | Structural equation model for the evaluation of national funding on R&D project of SMEs in consideration with MBNQA criteria[END_REF]. Under these circumstances, firms developing new technological skills may harvest the returns in new processes or new products, which lead to a competitive advantage [START_REF] Buisseret | What difference does it make? Additionality in the public support of R&D in large firms[END_REF][START_REF] Link | An econometric analysis of trends in research joint venture activity[END_REF].

Studies on the effect of R&D funding analyze the originality and profitability of the resulting innovations (Ebersberger, 2005;[START_REF] Guellec | The impact of public R&D expenditure on business R&D[END_REF][START_REF] Moore | R&D advancement, technology diffusion, and impact on evaluation of public R&D[END_REF]; the differential effects when using diverse innovation outcomes [START_REF] Mairesse | The importance of R&D for innovation: A reassessment using french survey data[END_REF]; and the mediating impact of industry and cooperation ( Ebersberger, 2005). These previous studies have as common finding the positive relationship between the R&D grants and the technological effect experienced by firms. Nonetheless, most of these studies focus their attention on the number of patents as an extended and consistent protection regime over countries, disregarding their potential shortcomings as R&D output measure. Firstly, patents serve as protection mechanisms that prevent imitation from competitors. Notwithstanding, they have a preemptive scope by reinforcing the first mover advantage, although many of these patents may never become actual products. Secondly, usually complex products require more than one patent, which makes them overestimate the innovation performance of the firm. Thirdly, patents grant protection to inventions lying on a limited definition of innovation according to the Oslo Manual; which leaves products and technologies of lesser inventive activity unprotected. Finally, patenting decisions usually entail strategic decisions, thus, firms may opt out from patents and retain their technology as trade secrecy. On the contrary, product innovations are the most tangible expressions of the techniques and inventive capacity of the firm. They reveal whether the company knows or dominates the technology, so as they embody the firms' knowledge in actual outcomes. They represent the creation of new or significantly improved ready-tomarket products, with new or improved technical characteristics or materials, or with improved functional characteristics, as well as products that encompass a simplified way of use (OECD, 2005). At last, firms leery of revealing their knowledge through patents, are not that modest when disclosing their new products. For these reasons and the pervasiveness of product innovations in SMEs in opposition to patents we claim that:

H.1: R&D grants beneficiary firms obtain more product innovations than that not receiving grant support.

CAPABILITIES AS INNOVATION DRIVERS

Public intervention supporting firm's innovation may encourage the development of new resources and capabilities within firms [START_REF] Buisseret | What difference does it make? Additionality in the public support of R&D in large firms[END_REF][START_REF] Link | An econometric analysis of trends in research joint venture activity[END_REF][START_REF] Torres-Barreto | Impacto tecnológico y económico de las ayudas públicas a la I+D+i: el papel de los recursos y las capacidades[END_REF], and our study also digs on them as firm characteristics that may influence the relationship between R&D funding and product innovations. Substantial bodies of literature claim that they must be unique to the firm to turn into the source of competitive advantage (Jay [START_REF] Barney | The resource-based view of the firm: Ten years after 1991[END_REF][START_REF] Helfat | The dynamic resource-based view: Capability lifecycles[END_REF][START_REF] Wernerfelt | A resource-based view of the firm[END_REF]. Therefore, what the firm knows, possess, does and integrates become critical elements in innovation. Mostly, resources and capabilities facilitate firms exploiting their inputs varyingly and more efficiently than its competitors, and transforming their knowledge into valuable products (Alvarez-Melgarejo & Torres-Barreto, 2018a, 2018c). Thus, resources and capabilities represent desirable static and dynamic skills when firms want to pursue new ventures with the auspice of public funds.

In particular, planning capabilities may provide firms with an explicit innovation strategy [START_REF] Hunter | Planning for innovation: A process oriented perspective[END_REF][START_REF] Petrick | Technology roadmapping in review: A tool for making sustainable new product development decisions[END_REF][START_REF] Phaal | Customizing the technology roadmapping approach[END_REF], and constitutes decision analysis resources of firms. On the risky road of innovation, to follow a determined innovation plan may turn out into a necessary guide to help firms to gather the benefits of R&D investments by minimizing risks [START_REF] Tomala | Innovation management: A synthesis of academic and industrial points of view[END_REF]; by motivating firms to move into yet unknown fields [START_REF] Clark | Leading product development[END_REF], driving opportunities of successful innovations. Also, previous research on key factors of invention´s success point out the necessity of a high level discipline of innovation planning [START_REF] Damanpour | The dynamics of the adoption of product and process innovations in organizations[END_REF], since rigorous tools and structured methods of planning innovation, together with the intervention of several areas within the company (marketing, engineering, design), are considered pillars of successful R&D outcomes, and competitive advantage producers [START_REF] Hunter | Planning for innovation: A process oriented perspective[END_REF]. Therefore, we posit that: H.2: R&D grant beneficiary firms obtain more product innovations if they follow an innovation plan.

Capabilities linked to human resources within firms have been considered to be a complementary asset to product innovation. From the resource based view perspective, an experienced human resource system creates and develops organizational capabilities that become the source of competitive advantages [START_REF] Saa-Perez | A resource-based view of human resource management and organizational capabilities development[END_REF]. Also, high experienced individuals are linked to creative success [START_REF] Chi | Selfexplanatirms: How students study and use examples in learning to solve problems[END_REF][START_REF] Hershey | The effects of expertise on financial problem solving: Evidence for goal-directed, problem-solving scripts[END_REF], and leadership and experience it's been proved to be linked to the number of ideas generated and the development of plans to implement them [START_REF] Vincent | Divergent Thinking, Intelligence, and Expertise: A Test of Alternative Models[END_REF].

But, although capabilities and experience of personnel are relevant in terms of obtaining innovation outcomes, the crux of the issue is the capability of firms to maintain a skilled and qualified human resource system that supports the innovation processes through an organizational culture focused on innovation [START_REF] Lado | Human resource systems and sustained competitive advantage: a competency-based perspective[END_REF][START_REF] Mcmahan | Alternative theoretical perspectives for strategic human resource management revisited: progress, problems and prospects[END_REF][START_REF] Muffatto | Corporate and individual competences: How do they match the innovation process[END_REF]. The capability of hiring and assuring the persistence of an R&D team attain competitive advantages for firms, whose organizational culture is focused on innovation, and have human resource systems with a higher propensity to create and develop new ideas [START_REF] Hunter | Planning for innovation: A process oriented perspective[END_REF][START_REF] Lau | The HR system, organizational culture, and product innovation[END_REF]. Therefore, we propose that: H.3: R&D grant beneficiary firms obtain more product innovations if they are capable of maintaining an updated R&D knowledge base.

Another key element considered in this research is the cooperation capability as mediator in the relationship between R&D funding and technological success. From the literature search, we found that a technological cooperation relationship may have positive as well as negative effects on firm´s innovation results. From the Organizational Learning Theory, cooperation among partners is understood as a mechanism that generates technological knowledge complementarities among partners; deriving into incremental and radical innovations. This is because a firm only possess some of the resources and capabilities needed in order to achieve technological success since those resources and capabilities are mostly tacit and complex and not easily obtained over time [START_REF] Pitelis | Edith Penrose's "The Theory of the Growth of the Firm" Fifty Years Later[END_REF], hence, making evident the importance of designing a proper R&D alliance portfolio [START_REF] Belderbos | Cooperative R&D and firm performance[END_REF][START_REF] Quintana | Relación tecnológica en los acuerdos de cooperación empresarial y generación de innovaciones[END_REF]. Cooperation may allow firm: to access supplementary technologies and assets, and also facilitates scale economies [START_REF] Baum | Don't Go It Alone: Alliance Network Composition and Startups' Performance in Canadian Biotechnology[END_REF][START_REF] Miotti | Co-operative R&D: why and with whom? An integrated framework of analysis[END_REF][START_REF] Mohnen | What type of enterprise forges close links with universities and government labs? Evidence from CIS 2[END_REF][START_REF] Teece | Towards an economic theory of the multiproduct firm[END_REF], to access to the partner's knowledge [START_REF] Belderbos | Cooperative R&D and firm performance[END_REF][START_REF] Fritsch | Innovation, regional knowledge spillovers and R&D cooperation[END_REF], to expand their knowledge capabilities [START_REF] Becker | R&D cooperation and innovation activities of firmsevidence for the German manufacturing industry[END_REF][START_REF] Kleinknecht | Why do firms cooperate on R&D? An empirical study[END_REF], to generate more tangible innovations when combining it with internal R&D efforts [START_REF] Cassiman | R&D cooperation and spillovers: some empirical evidence from Belgium[END_REF] and to increase the probability of getting new products as the number of cooperating partners raise [START_REF] Becker | R&D cooperation and innovation activities of firmsevidence for the German manufacturing industry[END_REF].

On the other hand, there is a number of authors postulating that competitive advantages linked to R&D cooperation have also disadvantages and risks associated. The empirical evidence in this sense is reduced and more studies in this field are believed to be necessary [START_REF] Ahuja | Collaboration networks, structural holes, and innovation: A longitudinal study[END_REF][START_REF] Powell | Interorganizational collaboration and the locus of innovation: Networks of learning in biotechnology[END_REF][START_REF] Stuart | Interorganizational alliances and the performance of firms: a study of growth and innovation rates in a high-technology industry[END_REF]. Opportunistic behaviors [START_REF] Das | Resource and risk management in the strategic alliance making process[END_REF], competition instead of cooperation [START_REF] Hamel | Competition for competence and interpartner learning within international strategic alliances[END_REF], extra costs of cooperation [START_REF] Doz | The evolution of cooperation in strategic alliances: Initial conditions or learning processes?[END_REF] and interdependence among cooperating partners [START_REF] Kale | Learning and protection of proprietary assets in strategic alliances: Building relational capital[END_REF] have been pointed out in the literature as main negative consequences of cooperation.

With different lines of thinking, we consider that studying the effects of cooperation on the relationship between R&D funding and innovative results of firms may be of a high interest, therefore we propose that: H.4: R&D grant beneficiary firms obtain more product innovations if they cooperate with partners for R&D activities.

Another capability considered in this research is the absorptive capability. This is the ability firms have to identify, assimilate and use knowledge from the environment (Alvarez-Melgarejo & Torres-Barreto, 2018b; [START_REF] Cohen | Innovation and Learning: The Two Faces of R & D[END_REF]. This includes not only the capability some firms have to copy and imitate innovations, but also the capability to identify basic research from third parties and turn it out into applied research by means of product innovations, utility models, or other forms of innovations that reach the market.

According to [START_REF] Cohen | Innovation and Learning: The Two Faces of R & D[END_REF], the absorptive capability represents an important part of firm´s skills needed to create new knowledge. Indeed, this capability is recognized to be a byproduct of previous investments firms have done regarding R&D, as the quantity of knowledge a firm accumulates through R&D processes, and the variety of this knowledge may influence the development of an absorptive capability by the side of firms. The more R&D investments a firm perform, the higher its capability to appreciate and value the external knowledge which is part of the absorptive capability, so, R&D investments definitely play an important role in the development of absorptive capability within firms by means of preparing them to assimilate technological knowledge, enabling firms to see different ways of perform activities or processes and finally, opening the door for the obtention of new innovations [START_REF] Cohen | Innovation and Learning: The Two Faces of R & D[END_REF].

The importance of absorptive capability has been explained by some authors in terms of the ability firms have to effectively transform the knowledge obtained from outside, and convert it on ideas, inventions, innovative services or products that could represent revenues for firms [START_REF] Zahra | Absorptive capacity: A review, reconceptualization, and extension[END_REF], making evident the significance of the acquired capability. Some other researchers have stated that acquisition and use of external information is closely related to the development of new products that overtake its competitors [START_REF] Cassiman | In search of complementarity in innovation strategy: Internal R&D and external knowledge acquisition[END_REF]. In spite of these statements, the relationship between absorptive capability and the development of new inventions have scarcely been explored (Atuahene-Gima, 1992; [START_REF] Kumar | Differential learning and interaction in alliance dynamics: A process and outcome discrepancy model[END_REF]. For this reason, we dig deeply on the precise ways a firm has to generate absorptive capability and found that a firm internally performing R&D activities is more prone to develop such a capability [START_REF] Mowery | Strategic alliances and interfirm knowledge transfer[END_REF][START_REF] Tilton | International diffusion of technology: The case of semiconductors[END_REF]. Also technological training activities contributes to boost the absorptive capability (Cohen & Levinthal, 1990) and, performing manufacturing activities play an important role as well, since such operations prepares companies to recognize and take advantage of new product's market information [START_REF] Abernathy | Patterns of industrial innovation[END_REF]. To empirically test this relationship, we propose that: H.5: R&D grant beneficiary firms obtain more product innovations if they have developed an absorptive capability.

After this argumentation, we expect that the impact of public support policies follows a different pattern on the production of product innovations, depending on the capabilities firms may develop. Figure 1 shows the complete model of this research.

METHODOLOGY

The empirical strategy consists of 1) testing whether these resources and capabilities, in combination with the aids, can lead to new products 2) testing whether these resources and capabilities alone lead a company to obtain product innovations, or increases the chances of getting them, but that they do not take to him to obtain helps. That is to say: it is not because of having those resources and capabilities that companies obtain aid.

Our study relies on a sample of 3,350 firms included in the Spanish Survey of Business Strategies (ESEE). We analyze manufacturing companies with more than ten employees through a panel data analysis including years 1998-2008. Since our dependent variable is discrete we chose a negative binomial (NB) model, as it has been suggested in previous studies [START_REF] Hilbe | Negative binomial regression[END_REF][START_REF] Mckenzie | Some simple models for discrete variate time series 1[END_REF]. We test both the fixed and random effects variants of NB specification. Finally, and consistently with previous studies suggesting that there might be differences between regional and national R&D policies (Afcha Chávez, 2011b, 2011a), we estimate different equations depending on the origin of the funding. Secondly, we address endogeneity by choosing a two-stage least squares specification, identifying those resources and capabilities which can lead to obtaining national or regional R&D funding, but that do not explain new product development. We compare both results.

MAIN FINDINGS

We study the influence of public R&D funds on product innovations, and we separately analyze the funds' origin on their effect. Moreover, we explore the role of resources and capabilities in obtaining product innovations through public R&D funds. In particular, we study the resources of own R&D spending and the industrialization level of the firm. Alternatively, we review the capabilities of innovation planning, the constant update of their knowledge base, the technological cooperation, and the absorptive capacity.

Our results indicate that public R&D funding from regional entities has a positive and direct effect on product innovations, but it is not the case for the national R&D programs. On the side of resources, the private investment in R&D leads to greater number of products. However, the industrialization level of the firm is not significant. On the side of capabilities, the innovation planning, the constant update of the knowledge base, and cooperation agreements with customers and suppliers positively influence the number of product innovations obtained. Nevertheless, we could not find a moderating effect of capabilities on the relationship between R&D funding, and product innovations, except for the case of national grants, highlighting a potential switch to more science oriented partners.

IMPLICATIONS FOR MANAGERS, PRACTICIONERS AND POLICY MAKERS

This research shed light on the effectiveness of R&D public funding and points out the necessity of reinforcing resources and capabilities that boost their economic and technologic position. The main findings regarding technologic impact give a hint to identify the technological impact derived from the public investments in private R&D activities; as well as to characterize the role that resources and capabilities of firms play in product innovation.

From our analysis, we can assert that public funding for R&D activities has an impact on the technological output of Spanish firms, and that resources and capabilities of firms do play a role in obtaining a higher technologic production (Torres-Barreto, 2017). Our findings contribute to the Resource-Based View of the firm in the sense that resources and capabilities do matter to explain heterogeneity in innovation. Moreover, the heterogeneity in outcomes is the result of strategically diverse set of resources and capabilities (Torres-Barreto, Martínez, Meza-Ariza, & Molina, 2016).

Our findings have an application from a regulatory perspective because they may contribute to the design of technology policy, and reducing market failure. Since the design of the particular funding programs, has a significant effect on the private expenditure on R&D, by complementing or substituting the private R&D expenditure. On the other hand, practitioners would also extract relevant conclusions from these results as our study highlights the necessary resources and capabilities needed for both apply to public grants and introducing new products.

Our results shed light on the literature about the effectiveness of R&D public funding by measuring the tangible outputs of the firm. However, we also call to the necessity of developing and reinforcing resources and capabilities to potentiate the effect of R&D grants on the number of product innovations. Therefore, our conclusions depict scenarios where policymakers can ascertain the spillovers of public R&D funding, as well as to anticipate the potential funding results depending on the resources and capabilities of the grantees. Moreover, managers can identify the relevant resources and capabilities to receive a grant, but also to extract its full potential.

Nevertheless, it raises another question, are these resources and capabilities ex-ante characteristics to R&D grant application, or are the funds that help to develop these resources and capabilities?

Figure 1 :

 1 Figure 1: Research model.

El escenario de financiación de