Latest revision as of 12:50, 15 June 2018

Abstract

In this part of paper we shall extend the formulation proposed by Babuška and co-workers for robustness patch test, in quality assessment of error estimators, to more general cases of patch locations especially in three dimensional problems. This is performed first by finding an asymptotic finite element solution at interior parts of a problem with assumed smooth exact solution and then adding a correction part to obtain the solution near a kinked boundary irrespective of other boundary conditions at far ends of the domain. It has been shown that the solution corresponding to the correction part may be obtained in an spectral form by assuming a suitable proportionality relation between the nodal values of a mesh with repeatable pattern of macro-patches.

Having found the asymptotic finite element solution, the performance of error estimators may be examined. Although in this paper we focus on the asymptotic behavior of error estimators, the method described in this part may be used to obtain finite element solution for two/three dimensional unbounded heat/elasticity problems with homogeneous differential equations. Some numerical results are presented to show the validity and performance of the proposed method.

Keywords: Recovery, Error Estimate, Superconvergent, Robustness test, three dimensional, Unbounded domains

1. INTRODUCTION

Efforts by pioneering engineers in late twentieth century led to understanding the necessity of error evaluation and use of adaptive procedures in finite element solutions. Nowadays, most of commercial finite element programs are equipped with tools for error estimation and mesh generation. Nevertheless, the need for further investigation is still felt by scientists [1].

The knowledge of error estimation in two dimensional problems has now reached to a level of maturity [2]. However, advances in engineering and technology have led to increase of demands for three dimensional analyses and thus it is of great interest to know whether the error analyses still work well in three dimensions. Although many of the concepts introduced in two dimensional problems may readily be applied to three dimensions, there are many especial cases that may not comprehensively be covered by the counterparts studied in two dimensions. Quality assessment of error estimators may be considered as one of the cases.

There are two basic approaches employed in studies for assessment of error estimators. The most classic approach, used almost since the first thoughts of FEM, is benchmarking. This category of approach was employed in many early studies on error estimation, see [3,4] for example, and it is still being used as a tool for giving some evidences of good performance, e.g. [5,6]. Although benchmark problems in two and three dimensional spaces have much in common, dealing with the three dimensional ones may need some more elaboration. For example behavior of error estimators near an edge in 3D problem, where two intersecting boundary exist, may not easily be understood from the results obtained at edges of a 2D benchmark. In fact an edge in 3D problem has much in common with a corner in a 2D problem, yet not quite similar.

The second approach is employed for studies on quality of approximate error in asymptotic manner, i.e. when the element sizes tend to zero, and has roots in basic assumptions of finite element method. Although this type of assessment approach was inherently used by some researches to give proofs for theorems and lemmas regarding the convergence of the proposed methods, e.g. see [7], for the first time a systematic computer based type of the approach was introduced by Babuška et al in a series of papers [8-10]. Other researchers followed some similar approaches to explore superconvergence in plate problems [11, 12].

The method was basically introduced for assessment of error estimators applied at interior parts of a fictitious unbounded domain and has sound mathematical bases. The basic feature of the method is using meshes with repeatable patterns of patches in order to reduce the size of the space needed for the finite element solution. The test procedure is performed on a solution found for a single patch and thus, in this respect, resembles to the patch tests usually used in element technology. Since the reduction of the size of the domain, i.e. from an infinite size to a patch with finite size, is permitted under some assumptions for asymptotic behavior of the exact solution, the test results are interpreted in an asymptotic point of view. The formulation given for two dimensional problems is general and may be applied to three dimensional problems with no much effort.

The method was further extended to patch locations near boundaries in two dimensional problems [13]. The new formulation includes a boundary layer solution, in an asymptotic manner, for a fictitious unbounded domain. However, the modified version of the formulation can only be applied to a domain with single flat boundary. This means that the method may not be applied for patches near corners in two dimensional cases. Of course, this is of not much importance since the frequency of dealing with such situations is very small compared to other cases studied in two dimensional cases.

For three dimensional problems, however, the modified version can only be used for patch locations near a flat boundary. This means that the cases in which the patch is located near two intersecting flat boundaries, i.e. an edge, cannot be studied. It is obvious that frequency of such cases is not comparable with that of corners in two dimensional problems. However, if the problem is solved for edges of a three dimensional body, the formulation may also be applied for corners of a two dimensional problem.

In this paper we shall extend the formulation to the situations that the patch is near an edge/corner of the domain. We shall employ a spectral formulation in which a proportionality property plays a key role. It will be shown that the proportionality property has roots in the nature of the governing differential equation. As will be shown later in this paper, the governing equation for boundary layer solution in integral form is of a homogeneous type. The formulation given in this paper is not only useful for finding asymptotic boundary layer solution but also applicable to many other unbounded problems with governing equations of similar type. This is an important feature of the method. Having found the formulation for the boundary layer solution near kinked boundaries, the results will be verified by a sample problem with exact solution on an unbounded domain. The test results and the discussion for quality of the error estimates are given in the next part of the paper.

The layout of the paper is as follows; In Section 2 we shall prepare some preliminary formulations and definitions to be used in subsequent sections. In Section 3 some basic assumptions initially made in [8-10] and basic relations for the test at interior parts will be revisited. The relations will be frequently referred to especially when we describe the step by step test procedure in its generalized form. Having recalled the basic relations, the formulations used for flat boundaries, introduced in [13], will be given in brief followed by our proposed method which is to be applied for kinked boundaries. In Section 4 we shall present a sample problem to demonstrate that the proposed method is able to produce finite element solution for general unbounded domain with tractions or prescribed values at boundaries exhibiting decay condition. We shall conclude our discussions in Section 5.

2. PRELIMINARIES

In this section we prepare the preliminary tools for study of asymptotic behaviors of finite element solutions as well as error estimation schemes.

The problems of interest in this study are elliptic ones, e.g. heat transfer or elasticity problems, with generic differential equation/equations as

(1-a)

to be solved on ${\textstyle \Omega }$ with following boundary condition

on

(1-b)

and

on

(1-c)

Where ${\textstyle {\Gamma }_{u}\cap {\Gamma }_{t}=\varnothing }$ and ${\textstyle {\Gamma }_{u}\cup {\Gamma }_{t}=\partial \Omega }$ . In above ${\textstyle {\boldsymbol {S}}}$ is a suitable operator for defining the gradients of ${\textstyle {\boldsymbol {u}}}$ , ${\textstyle {\boldsymbol {D}}}$ is a suitable matrix for modulus of material and ${\textstyle {\boldsymbol {t}}}$ is a traction vector associated with the problem type. Also ${\textstyle {\boldsymbol {\tilde {n}}}}$ is a matrix containing components of unit normal to the boundary arranged in an appropriate form for heat or elasticity problems.

Discretization of the domain to a series of elements and approximation of ${\textstyle {\boldsymbol {u}}}$ as ${\textstyle {\boldsymbol {u}}\cong {\boldsymbol {u}}_{h}={\boldsymbol {N}}{\boldsymbol {\overline {u}}}_{h}}$ and application of Galerkin type of weighted residual formulation in weak form (or equivalent minimization of a suitable functional) leads to following well known matrix equation

(2)

Where ${\textstyle {\boldsymbol {N}}}$ is the set of shape functions used and also ${\textstyle {\boldsymbol {B}}\equiv {\boldsymbol {SN}}}$ . It is noteworthy to remember that the following relation always exists between the exact solution and the finite element one

(3)

Generally there will be a discrepancy between the exact solution, which is not available for most problems, and the finite element solution. This deviation constitutes the error at each point

(4-a)

The error distribution can also be expressed in terms of the gradients

or

(4-b)

Where ${\textstyle {\sigma }_{ex}={\boldsymbol {D}}{\epsilon }_{ex}={\boldsymbol {DS}}{\boldsymbol {u}}_{ex}}$ and ${\textstyle {\sigma }_{h}={\boldsymbol {D}}{\epsilon }_{h}={\boldsymbol {DB}}{\boldsymbol {\overline {u}}}_{h}}$ .

If the shape functions are chosen properly, reduction of the maximum mesh size must result to reduction of errors, i.e. ${\textstyle {\boldsymbol {e}}\rightarrow {\boldsymbol {0}}}$ as ${\textstyle h_{max}\rightarrow 0}$ . This is usually referred to as completeness of the set of shape functions.

It is a priori understood that in absence of any singularity in the domain, the ideal rate for convergence of the finite element errors is of order ${\textstyle O(h^{P+1})}$ for ${\textstyle {\boldsymbol {e}}_{u}}$ and ${\textstyle O(h^{P})}$ for ${\textstyle {\boldsymbol {e}}_{\epsilon }}$ and ${\textstyle {\boldsymbol {e}}_{\sigma }}$ where ${\textstyle h}$ is an average mesh size for the solution. This simple information constitutes a basis for study of error behavior and evaluation of error estimates in asymptotic form.

2.1 Error estimation

Approximate evaluation of error can be performed through estimation of its absolute value, irrespective of the signs, or calculation of the difference between a set of enhanced answers and those of finite element method. Most of residual based error estimators fall within the former approach and almost all recovery based error estimators employ the latter approach.

In an error analysis a suitable norm is usually defined to give an indication for the exact error level over the whole domain, e.g. ${\textstyle L_{2}}$ norm of errors

or

(5)

Here ${\textstyle {\Omega }^{'}}$ represents a partition of the domain where the error analysis is to be performed. When evaluation of an error estimator is required, benchmark problems with exact solution are selected and the true error and the approximate one are compared. Having evaluated approximate errors, i.e. ${\textstyle {\boldsymbol {e}}_{u}^{\ast }}$ , ${\textstyle {\boldsymbol {e}}_{\epsilon }^{\ast }}$ or ${\textstyle {\boldsymbol {e}}_{\sigma }^{\ast }}$ , a dimensionless parameter as the ratio of the error values, called “effectivity index”, is defined to show the quality of error estimates

(6)

where

(7)

Obviously when ${\textstyle \theta }$ is close to unity, the performance of the estimator is considered satisfactory. It is expected that the effectivity index tend to unity when the mesh size is reduced sufficiently, i.e. ${\textstyle \theta \rightarrow 1}$ as ${\textstyle h\rightarrow 0}$ . If such a property exists, the error estimator is called “asymptotically exact” [14]. It has been shown that only few estimators exhibit such a property (the reader may consult to references 3,4 and 10).

However, in study of asymptotic behavior of error estimators, it happens that the asymptotic effectivity indices take values close to unity. When comparison of different types of estimators is aimed, an error estimator is said to be more robust if the effectivity index calculated is closer to unity.

3. ROBUSTNESS PATCH TEST

In this section we describe the test introduced by Babuška et al in 1994. Although the procedure can be found, with many proofs for theorems and lemmas, in references [8] and [9] we repeat the core of the formulation here since it is needed for generalization of the formulation to wider set of cases.

The test is basically a patch test since all computations at last is performed on a small group of elements with specific boundary conditions. The whole procedure resembles to the patch tests usually employed for examining completeness of shape functions in element technology. The reader can find the basis and history of such patch tests in reference [15]. But before explaining the robustness test, it is worthy to recall the basic features of the patch test in element technology in order to comprehend similarities and differences with the robustness patch test for error estimators.

It is well understood that if in a simple patch test an element can reproduce an exact field with order similar to the order of the shape functions then the asymptotic rate of convergence of the solution is guaranteed. The test is in fact a finite element solution on a group of elements with boundary conditions extracted from the exact field. The choice of boundary condition for such a test is trivial to determine and this is usually accomplished by considering prescribed values for main function, i.e. displacements, or its derivatives as prescribed tractions.

It may be thought that similar strategy is applicable when asymptotic convergence of an error estimate, e.g. a superconvergent one, is to be studied. It is true in some respects but the test will not be a complete one. It is possible to construct a test on a group of elements with prescribed boundary conditions extracted from a polynomial, as an exact solution, with one order higher than that of the shape functions. In such a study the aim will be evaluation of capability of reproducing the exact errors. But we note that the boundary conditions considered are not necessarily similar to those which exist between the patches of elements in a real finite element solution. Therefore the results obtained from such a patch test may be considered just as a rough indication for capability of the error estimator.

The deficiency is removed in the patch test proposed by Babuška et al in an elegant manner. The test is constructed by some basic assumptions for mesh configuration and smoothness of the solution. The robustness of an error estimation, in a norm as (5) for instance, is defined by evaluation of upper and lower bounds, ${\textstyle {\theta }_{L}}$ and ${\textstyle {\theta }_{U}}$ , for effectivity index of (6). Therefore the following index is defined

(8)

If the error estimator precisely reproduces the error then ${\textstyle {\theta }_{L}={\theta }_{U}=1}$ and thus ${\textstyle R=0}$ . However for practical use, robustness values between ${\textstyle 0.8<R<1.2}$ are considered ideal. Since the assessment is performed in an asymptotic manner, the exact solution for evaluation of the exact error is considered as a polynomial of one order higher than that of shape functions representing largest terms in a Taylor expansion of error when the size of elements, at limit, is infinitesimal.

3.1 Assumptions

The fundamental assumptions are as follows

It is assumed that in all directions the mesh is locally uniform and is composed of a series of repeatable pasterns, as the basic cell with dimension ${\textstyle h}$ especially inside ${\textstyle {\Omega }_{1}}$ in Figure 1. A single cell of the repeatable pattern is considered as a patch for the test. The reader may notice that such patches have different definitions from the ones used for recovery methods. Note that no assumption is made about the mesh outside of ${\textstyle {\Omega }_{1}}$ .


(a)	(b)

Figure 1. The domain, subdomain and basic cell considered for robustness patch test; (a) Three dimensional presentation of the domain; (b) Projection of the domain on a two dimensional plane.

Study on asymptotic convergence of the solution requires that the elements sizes be close to zero and so the dimensions of patches consisting a number of such elements like patches shown in Figure 1, i.e. basic cell with dimension ${\textstyle h}$ and the larger patches with dimensions ${\textstyle h_{0}}$ and ${\textstyle h_{1}}$ ( ${\textstyle h_{0}}$ and ${\textstyle h_{1}}$ are multiples of ${\textstyle h}$ ). In this regard it is also assumed that that ${\textstyle h}$ tends to zero with a faster rate than ${\textstyle h_{0}}$ (and ${\textstyle h_{1}}$ ). Moreover, it is assumed that the basic cells are cuboid (in three dimensions) with edges parallel to axes, though this is not necessary in general.

It is assumed that the exact solution is smooth at vicinity of the patch, i.e. inside macro patch ${\textstyle {\Omega }_{1}}$ . This means that the exact solution has bounded ${\textstyle \left(p+2\right)}$ derivatives at the location of the patch (or in other words the patch is far from singular points). Besides, not all of the derivatives with order of ${\textstyle \left(p+1\right)}$ vanish at the location of the patch.

The mesh is sufficiently refined outside ${\textstyle {\Omega }_{1}}$ such that the error converges almost with the optimal rate in ${\textstyle L_{2}}$ norm in ${\textstyle {\Omega }_{1}}$ , i.e.:

${\textstyle {\Vert {\boldsymbol {e}}_{u}\Vert }_{L_{2},{\Omega }_{1}}\leq Ch^{p+1-\lambda }}$

where ${\textstyle C}$ is a constant independent of ${\textstyle h}$ , and ${\textstyle \lambda }$ representing the power of singularity, is sufficiently small. This means that the pollution error is negligible. This assumption practically holds for meshes created by adaptive algorithms which try to distribute, equally, the energy norm of the error over the elements.

3.2 Basic features of the test

Before focusing on the error, asymptotic finite element solutions should be obtained. This means that the problem of (1) should be solved either irrespective of the main boundary conditions, when the patch is located inside, or just with condition specified when the patch located close to a boundary but still irrespective of other boundary conditions. The information available is that the exact solution, in asymptotic sense, is of a polynomial form. If such a finite element solution is obtained the results are generalized to many problems with similar type and different conditions at boundaries. Error evaluation comes after determining the finite element solution.

From the smoothness assumptions it is concluded that the exact solution may be written in the form of a Taylor series around a point inside the macro-patch ${\textstyle {\Omega }_{0}}$ as

(9)

with ${\textstyle a_{i}}$ being appropriate coefficients and ${\textstyle \Delta x=x-x_{0}}$ and so forth for ${\textstyle \Delta y}$ and ${\textstyle \Delta z}$ where ${\textstyle x_{0}}$ , ${\textstyle y_{0}}$ and ${\textstyle z_{0}}$ are coordinates of a point inside the patch. With the knowledge of asymptotic convergence of the finite element solution it follows that it is just necessary to consider a complete polynomial of order ${\textstyle \left(p+1\right)}$ as the asymptotic exact solution. For scalar problems, for instance

,

.

(10)

Here ${\textstyle X}$ , ${\textstyle Y}$ and ${\textstyle Z}$ are suitable normalized coordinates. For testing the estimates on linear elements in three dimensional elasticity problems, for example, we write the polynomial vector as

,

(11)

where ${\textstyle {\boldsymbol {1}}_{1}={\langle {\begin{array}{ccc}1&0&0\end{array}}\rangle }^{T}}$ , ${\textstyle {\boldsymbol {1}}_{2}={\langle {\begin{array}{ccc}0&1&0\end{array}}\rangle }^{T}}$ and ${\textstyle {\boldsymbol {1}}_{3}={\langle {\begin{array}{ccc}0&0&1\end{array}}\rangle }^{T}}$ .

Now with the aid of the following trivial relation

(12)

in which ${\textstyle {\boldsymbol {\overline {u}}}_{ex}}$ is a vector of exact values at nodal points of the mesh, one may use Equation (3) and write

(13)

Wherever the mesh is repeated, the right hand side of (13) is a periodic function since ${\textstyle {\boldsymbol {u}}_{ex}-{\boldsymbol {N}}{\boldsymbol {\overline {u}}}_{ex}}$ is periodic. The proof is very easy and can be accomplished by writing the monomials of (10) or (11) locally in each element and substituting in expression of ${\textstyle {\boldsymbol {u}}_{ex}-{\boldsymbol {N}}{\boldsymbol {\overline {u}}}_{ex}}$ . It follows that the new function is not dependent on the position of the origin of the local coordinates and therefore if the element pattern is repeated the function will be repeated as well (see [16] for more explanations).

The repeatability property of the right hand side of Equation (13) plays a key role in reducing the size of the problem to be solved. The periodic nature of the right hand side of (13) leads to repeatability of the left hand side of the equation. Therefore one may write the equation on a smallest repeatable unit of the mesh, or on a larger repeatable portion consisting of the units, here known as ${\textstyle {\Omega }^{'}}$

(14)

in which ${\boldsymbol {\overline {u}}}_{h}^{p}={\boldsymbol {\overline {u}}}_{h}-$ ${\boldsymbol {\overline {u}}}_{ex}$ . Equation (14) should of course be solved with a set of periodic boundary conditions. In the next subsection we shall see that depending on the patch location ${\textstyle {\Omega }^{'}}$ may be interpreted as a patch volume itself or volume of a column of patches normal to the boundary.

It can be seen that Equation (14) is insensitive to constant values for ${\textstyle {\boldsymbol {\overline {u}}}_{h}^{p}}$ , i.e. rigid body motion of the solution. This implies that after finding ${\textstyle {\boldsymbol {\overline {u}}}_{h}^{p}}$ , as will be described later, the asymptotic finite element solution is obtained as

(15)

Where ${\textstyle {\boldsymbol {C}}}$ is a vector containing constant values and is evaluated from equivalence of mean values for the exact and finite element solutions [9, 10, 13, 16]

(16)

Having found the asymptotic finite element solution from ${\textstyle {\boldsymbol {u}}_{h}={\boldsymbol {N}}{\boldsymbol {\overline {u}}}_{h}}$ , the exact norm of errors and that of estimated ones are calculated. Since the exact solution is expressed in terms of monomials with unknown coefficients ${\textstyle {\boldsymbol {a}}}$ , as Equation (10) or (11), the norms of errors in squared form can be written as

(17)

With analogy, the approximate error may also be written as

(18)

This makes it easy to find the upper and lower bounds of the squared effectivity index

(19)

by evaluation of eigenvalues of the following matrix

(20)

Where the matrix ${\textstyle {\boldsymbol {H}}}$ is defined as

,

(21)

In above equation ${\textstyle \Psi }$ and ${\textstyle \Lambda }$ are normalized eigenvectors and diagonal matrix of eigenvalues of ${\textstyle {\boldsymbol {E}}_{0}}$ , respectively. Hence

,

(22)

It should be noted that relations (12) to (16) are useful for reducing the size of the problem only when a sort of periodicity exists in the solution. This will be the case as long as the exact solution is expressed in terms of monomials of order ${\textstyle p+1}$ and the mesh pattern is repeatable. When the patch is far enough from the main boundary conditions, i.e. located at interior of the domain, the periodicity may be assumed for all directions. This allows us to reduce the size of the problem to a domain as small as a patch. In the case that patch is located near a boundary, the repeatability of the mesh is lost in normal direction to the boundary and thus the solution loses its periodicity in that direction. This latter case needs additional treatment which we shall briefly explain in the forthcoming subsections. In the case that the patch is located at a corner of the main domain, the periodicity will be lost in all directions. We shall also refer to this case after revisiting the patch test near a single flat boundary.

3.3 Patch test at interior parts of the domain

As mentioned above, when patch is located at interior parts of the domain, the patch pattern is considered repeatable in all directions. Now since the source term in the right hand side of (14) is repeatable, the finite element solution will be periodic. This allows us to write Equation (14) just on a single patch as

(23)

which is to be solved with periodic boundary conditions i.e.

${\boldsymbol {\overline {u}}}_{h(Z+L_{z})}^{p}={\boldsymbol {\overline {u}}}_{h(Z)}^{p}$

(24)

with ${\textstyle L_{x}}$ , ${\textstyle L_{y}}$ and ${\textstyle L_{z}}$ being the patch dimensions (period length) along the three directions. Having found the periodic solutions, the asymptotic finite element solution is evaluated from (15) and (16). The test procedure will be completed by estimating the errors and calculation of upper and lower bounds of the effectivity index as indicated in Equations (17) to (22).

3.4 Patch test at parts near flat boundaries

To extend the formulation to more general cases, here we revisit the procedure near flat boundaries. Now supposing that a flat boundary exists at ${\textstyle z=0}$ and the patch pattern is repeatable in all directions. Equation (14) may be rewritten as

(25)

Clearly the periodicity of solution is lost along z axis since the repeatability of the patch pattern is terminated at ${\textstyle z=0}$ although the source term is still periodic. Hence an additional solution is considered, as boundary layer type, to those obtained for interior parts

(26)

Substituting (26) in (25) and noting that (14) holds for solutions pertaining to interior patches, then provided that no extra forces is imposed from boundaries, i.e. when dealing with Dirichlet boundary conditions, the following integral equation is obtained

(27)

It is clear that as we go far from the boundary, i.e. ${\textstyle z\rightarrow \infty }$ , we should obtain ${\textstyle {\boldsymbol {\overline {u}}}_{hB}^{p}\rightarrow {\boldsymbol {\overline {u}}}_{h}^{p}+}$ ${\boldsymbol {C}}$ . Therefore Equation (27) must be solved with two periodic conditions along ${\textstyle x}$ and ${\textstyle y}$ and following decay condition

(or

) when

(28)

Where ${\textstyle d_{z}}$ denotes the distance along ${\textstyle z}$ . It may be noted that decay condition on ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ is stronger than the one defined for its gradients, i.e. ${\textstyle {\boldsymbol {B}}{\boldsymbol {\overline {u}}}_{h0}^{p}}$ . In fact application of decay condition in latter form may lead to constant values of ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ for large values of ${\textstyle z}$ which is still consistent with formulation for interior parts since constant values are considered a priori. This is true basically when we are dealing with asymptotic finite element solutions in which the solutions for two fictitious domains are added together as Equation (26). However, when a general unbounded domain is to be solved, constant values for the solution at far ends of the domain may have no physical meaning and thus in such situations the first form of decay condition must be satisfied rigorously. In a sample problem we shall give solution to an unbounded domain in order to verify the method to be described in sequel.

The boundary conditions considered in this study are

1. Dirichlet boundary condition: In this type we apply

(29)

Substitution of (26) instead of the last two terms in right hand of (15) as

(30)

and application of (29) leads to following condition

(31)

From which it follows that

(32)

If homogeneous Dirichlet condition is of concern then the monomials in (10) or (11) should be

chosen accordingly.

2. Neumann boundary condition: In this type of condition the tractions extracted from the exact solution are applied to the boundary of the domain, i.e.

(33)

to be substituted in (1-c). If Equation (2) is rewritten with ${\textstyle {\boldsymbol {\overline {u}}}_{h}}$ replaced with its equivalent from (26) and (27) then

(34)

Where

(35)

It follows that

(36)

In fact we only need to consider elements near boundary since Equation (27) holds for interior parts. Therefore

(37)

In above ${\textstyle {{\Omega }^{'}}_{bnd}}$ denotes volume of all elements in the column having at least one node on the boundary, ${\textstyle {\left({\boldsymbol {\overline {u}}}_{h0}^{p}\right)}_{bnd}}$ denotes associative subset of ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ , and ${\textstyle {\boldsymbol {F}}_{prs}}$ represents the right hand side of Equation (36).

Now we are ready to apply the boundary condition to a problem asymptotically modeled with a periodic mesh pattern and a single flat boundary. Considering a 3D mesh of patches like Figure 2-a, we assume that the boundary exists at plane ${\textstyle z=0}$ . The mesh pattern is considered repeatable along ${\textstyle x}$ , ${\textstyle y}$ and ${\textstyle z}$ , with ${\textstyle -\infty <x<\infty }$ , ${\textstyle -\infty <y<\infty }$ and ${\textstyle 0<z<\infty }$ . A column of patches along ${\textstyle z}$ is selected. For simplicity we use Figure 2-b instead of Figure 2-a as a two dimensional presentation of the three dimensional column. Note that we are aiming at solving Equation (27) or (37) in which ${\textstyle {\Omega }^{'}}$ is volume of the column.


(a)	(b)

Figure 2. A column with infinite length; (a) Three dimensional presentation of the column with one end at boundary consisting of repeatable patterns as basic cells, (b) Two dimensional presentation of the infinite column with repeatable pattern and the numbering used for different levels.

After numbering the patch levels consecutively towards the middle of the domain as shown in Figure 2-b, in the matrix equation of each patch the degrees of freedoms not falling on interfaces of patches are statically condensed. The matrix equation of the asymptotic finite element solution pertaining to the remaining degrees of freedoms is assembled as

{\left[{\begin{array}{ccccccccc}{\boldsymbol {k}}_{0,0}&{\boldsymbol {k}}_{0,1}&&&&&&&\\{\boldsymbol {k}}_{1,0}&{\boldsymbol {k}}_{1,1}&{\boldsymbol {k}}_{1,2}&&&&&&\\&{\boldsymbol {k}}_{2,1}&{\boldsymbol {k}}_{2,2}&{\boldsymbol {k}}_{2,3}&&&&&\\&&&&\ddots &&&&\\&&&&&{\boldsymbol {k}}_{n,n-1}&{\boldsymbol {k}}_{n,n}&{\boldsymbol {k}}_{n,n+1}&\\&&&&&&&&\ddots \end{array}}\right]}_{\infty \times \infty }\cdot {\left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,0}\\{\boldsymbol {\overline {u}}}_{h0}^{p,1}\\{\boldsymbol {\overline {u}}}_{h0}^{p,2}\\\vdots \\{\boldsymbol {\overline {u}}}_{h0}^{p,n-1}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n+1}\\\vdots \end{array}}\right\}}_{\infty \times 1}={\left\{{\begin{array}{c}{\boldsymbol {F}}_{0}\\{\boldsymbol {0}}\\{\boldsymbol {0}}\\\vdots \\{\boldsymbol {0}}\\{\boldsymbol {0}}\\{\boldsymbol {0}}\\\vdots \end{array}}\right\}}_{\infty \times 1}

(38)

Where ${\textstyle {\boldsymbol {k}}_{m,n}={\boldsymbol {k}}_{n,m}^{T}}$ is the relative stiffness matrix of ${\textstyle m}$ th and ${\textstyle n}$ th levels, and ${\boldsymbol {\overline {u}}}_{h0}^{p,k}$ is the set of nodal values at level ${\textstyle k}$ .

It can be seen that each row of the matrix represents the relation between levels ${\textstyle n-1}$ , ${\textstyle n}$ and ${\textstyle n+1}$ except for the first level. Also, since the patch pattern is the same for all levels, the out of diagonal matrices, below and above the diagonal, are similar. Extracting the ${\textstyle n}$ th row of the matrix Equation (38), one can show it as

(39)

By appropriate arrangement we have

\left[{\begin{array}{cc}{\boldsymbol {A}}&{\boldsymbol {B}}\\{\boldsymbol {0}}&{\boldsymbol {I}}\end{array}}\right]\cdot \left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,n+1}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n}\end{array}}\right\}=\left[{\begin{array}{cc}{\boldsymbol {0}}&-{\boldsymbol {C}}\\{\boldsymbol {I}}&{\boldsymbol {0}}\end{array}}\right]\cdot \left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,n}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n-1}\end{array}}\right\}

(40)

Supposing that the coefficient matrix in the left hand side of (40) is invertible, it follows that

\left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,n+1}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n}\end{array}}\right\}=\left[{\begin{array}{cc}-{\boldsymbol {A}}^{-1}{\boldsymbol {B}}&-{\boldsymbol {A}}^{-1}{\boldsymbol {C}}\\{\boldsymbol {I}}&{\boldsymbol {0}}\end{array}}\right]\cdot \left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,n}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n-1}\end{array}}\right\}

(41)

One may write a relation between two sets of such defined vectors for two levels, say levels ${\textstyle m}$ and ${\textstyle n}$ , as

\left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,n+1}\\{\boldsymbol {\overline {u}}}_{h0}^{p,n}\end{array}}\right\}={\left[{\begin{array}{cc}-{\boldsymbol {A}}^{-1}{\boldsymbol {B}}&-{\boldsymbol {A}}^{-1}{\boldsymbol {C}}\\{\boldsymbol {I}}&{\boldsymbol {0}}\end{array}}\right]}^{n-m}\cdot \left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{h0}^{p,m}\\{\boldsymbol {\overline {u}}}_{h0}^{p,m-1}\end{array}}\right\}

(42)

Solution of above recursive relation can be obtained by transformation of the variables to a space with bases of normalized eigenvectors of the coefficient matrix in Equation (41). Invariably, the eigenvalue problem may be solved using matrices in Equation (40), as a generalized eigenvalue problem, to avoid problems encountering for inversion of ${\textstyle {\boldsymbol {A}}}$ . Therefore following generalized eigenequation is defined

or

(43)

Solution of the eigenequation of (43) results in a series of eigenvalues ${\textstyle {\mu }_{i}}$ and eigenvectors ${\textstyle {\boldsymbol {w}}_{i}}$ , as many as two times of DOFs at each level, from which those with decay property, as ${\textstyle \vert {\mu }_{i}\vert \leq 1}$ , are nominated for construction of the asymptotic finite element solution. We note that in the new space constructed by series of ${\textstyle {\boldsymbol {w}}_{i}}$ the recursive relation for each base vector at two different levels, ${\textstyle m}$ and ${\textstyle n}$ , is written as

(44)

Also each eigenvector is of the form

(45)

Now assume that the column with infinite length is cut and the boundary conditions are to be imposed at current level with nodal values as ${\textstyle {\varphi }_{i}}$ . If the levels are renumbered accordingly then for ${\textstyle n}$ levels farther the correction to the periodic finite element solution is evaluated as

(46)

In which ${\textstyle n_{c}}$ is the number of DOFs at each level and ${\textstyle c_{i}}$ denotes the unknown coefficient for the mode number ${\textstyle i}$ to be determined from the boundary condition specified in (31) or (37).

Having found the correction required for boundaries, ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ , the periodic finite element solution may be calculated as Equation (26). Finally, the asymptotic finite element solution is obtained as

(47)

Now one can apply the error estimation procedure to study its asymptotic convergence behavior.

As is seen the formulation is suitable for flat boundaries only. In case that the patch is located near a kinked boundary the periodicity of the solution is lost along to directions normal to the edge line. The problem becomes worse when a corner patch near intersection of three planes, is considered. In this latter case the periodicity is lost along all directions. In the next section we shall extend the formulation to such wider range of applications.

3.5 Generalization of the test for application to patches near kinked boundaries

In this section we shall give a more general formulation for evaluation of the asymptotic finite element solutions near edges or corners of a three dimensional domain.

We assume that two intersecting planes exist at ${\textstyle x=0}$ and ${\textstyle y=0}$ making a straight edge along ${\textstyle z}$ axis (Figure 3). A repeatable mesh pattern is again constructed. Since the problem is studied in asymptotic form, the mesh is considered to be consisting of infinite number of patches along the three axes. Equation (25) is again written

(48)


(a)	(b)

Figure 3. A typical case for patch location near two intersecting flat boundaries; (a) The domain shaded in gray is an unbounded domain along ${\textstyle x}$ and ${\textstyle y}$ axes and is of a single-cell-width thickness, (b) Two dimensional presentation of patch location near kinked boundaries which is similar to a corner point.

Here ${\textstyle {\Omega }^{'}}$ represents an infinite domain with thickness of a patch along ${\textstyle z}$ axes since the periodicity of the solution retains just along ${\textstyle z}$ . To avoid confronting difficulties for three dimensional presentation of the problem, we look in to the problem as a two dimensional one, i.e. its projection on ${\textstyle x}$ ${\textstyle y}$ plane (see Figure 3-b). Again we assume that the asymptotic periodic solution is consisting of two parts as Equation (26) substitution of which in (48) gives similar equation as (27). Clearly the solution must be found satisfying Dirichlet or Neumann boundary conditions similar to (30) or (37) and two decay conditions along ${\textstyle x}$ and ${\textstyle y}$ . Thus the statement of the problem is to solve

(49)

with boundary conditions (also see (32) and (37))

or

(50)

at ${\textstyle x=0}$ or ${\textstyle y=0}$ , periodic conditions along ${\textstyle z}$ and following decay conditions

(or

) when

and

(51)

Where here again ${\textstyle d_{x}}$ and ${\textstyle d_{y}}$ denote distances along ${\textstyle x}$ and ${\textstyle y}$ .

The key point of the solution is isolation of a patch from , as before, and writing suitable relations between the DOFs. Supposing that the solution may be found spectrally in an appropriate space, some proportionality relations are assumed between adjacent patches. The reason for assuming such relations may be explored in Appendix. For instance in Figure 4-a following relations are considered

,

(52-a)

for direction and

,

(52-b)

for ${\textstyle y}$ direction. Similar relation is valid for any two nodes with distance ${\textstyle L_{x}}$ or ${\textstyle L_{y}}$ for example

,

(52-c)

and so forth. In above ${\textstyle {\boldsymbol {\overline {u}}}}$ denotes spectral nodal values in a new space to be defined later.


(a)	(b)

Figure 4. A macro patch consisting of ${\textstyle 3\times 3}$ patches; (a) The proportionality relations are written as (52), between nodal values along both directions, (b) Schematic presentation of nodes in the macro patch nominated for satisfaction of the boundary condition.

Now the solution starts with writing all relations similar to (52) for a patch and all surrounding ones, i.e. nine patches in Figure 4, in a matrix form. It can be seen that all nodal values may be written in terms of few of them in central patch. Note that the nodal values inside the central patch must be condensed before writing the proportionality relations. Thus the remaining nodal values will be those located at two edges of patches, i.e. left and bottom edges. Hence if we focus on all similar DOFs at the two similar sides of the patches, the following matrix relation containing proportionality properties, like relations (52), can be written

(53)

Where ${\textstyle {\boldsymbol {\overline {u}}}_{D}}$ and ${\textstyle {\boldsymbol {\overline {u}}}_{I}}$ denote dependent and independent nodal values in macro patch of ${\textstyle 3\times 3}$ in Figure 4-a. The independent nodal values in Figure 4-a may be considered as those falling in the dash line area. The reader may notice that the relation may also be written in a macro patch. Here again the matrix equation obtained from (49) may be presented as (38). To avoid difficulties pertaining to numbering of the patches along two or three axes, we use an abstract notation with rearranged form as the following one

\left[{\begin{array}{cc}{\boldsymbol {K}}_{II}&{\boldsymbol {K}}_{ID}\\{\boldsymbol {K}}_{DI}&{\boldsymbol {K}}_{DD}\end{array}}\right]\cdot \left\{{\begin{array}{c}{\boldsymbol {\overline {u}}}_{I}\\{\boldsymbol {\overline {u}}}_{D}\end{array}}\right\}=\left\{{\begin{matrix}{\boldsymbol {0}}\\{\boldsymbol {f}}\end{matrix}}\right\}

(54)

The first row of the above matrix equation represents a generic row of a matrix equation similar to (38) and the second row represents the remaining of the equations most of which have zero right hand side, i.e. components of ${\textstyle {\boldsymbol {f}}}$ , except for the ones corresponding to nodes at Neumann boundaries. Substituting (53) in (54), the first set of equations gives

(55)

In which

(56)

Equation (55) means that ${\textstyle {\boldsymbol {\overline {u}}}_{I}}$ may be constructed from vectors in null space of ${\textstyle {\boldsymbol {Q^{'}}}}$ . To evaluate such null space the determinant of ${\textstyle {\boldsymbol {Q}}^{'}}$ is set to be zero

(57)

It may be observed that two unknowns, ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ , appear in one equation. This means that numerous pairs may be found satisfying Equation (57). To solve the problem one may choose a value for one unknown and find the other through the characteristic equation. Hence

(58)

To satisfy the decay condition (51) we define another function as

(59)

Generally, several relations like (59) exist since characteristic equation (57) possesses several roots for a given ${\textstyle {\mu }_{1}}$ . Therefore

when

,

(60)

Now supposing that ${\textstyle {\phi }_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{1})}}$ belongs to the null space of ${\textstyle {\boldsymbol {Q^{'}}}}$ , one may write the non trivial solution of (55) as

(61)

Where integration is taken over all feasible values of ${\textstyle {\mu }_{1}}$ . In Equation (61) coefficients ${\textstyle c_{i({\mu }_{1})}}$ are to be determined from boundary conditions at ${\textstyle x=0}$ or ${\textstyle y=0}$ . The number of inner summation is dependent on the number of ${\textstyle {\phi }_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{1})}}$ for which ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ satisfy (59) and may vary when ${\textstyle {\mu }_{1}}$ falls within different regions.

Since ${\textstyle {\mu }_{1}}$ may take complex values the symbolic integral of (61) is taken over a feasible area on Gaussian plane. Because the decay condition imposes the condition of ${\textstyle \vert {\mu }_{1}\vert \leq 1}$ , the feasible area is a circle with radios of unity and center at origin. However, we note that for problems with real differential equation, the exact solution and the finite element one should be real. This means that from values in feasible domain those leading to real answers should be selected implying that the complex values, and the associated vectors, should be in the form of conjugate pairs. Thus the integration in (61) should in fact be taken over real values for ${\textstyle {\mu }_{1}}$ while letting ${\textstyle {\mu }_{2}}$ to be complex. To ensure that all real values are covered by ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ Equation (61) may be rewritten as

(62)

or

(63)

where of course it is needed to find ${\textstyle {\mu }_{1}}$ in terms of ${\textstyle {\mu }_{2}}$ as Equation (58) and (59). Note that the negative values have been excluded from the integration domain because of two reasons. Firstly, change of the sign of the proportionality factor represents an oscillatory solution, i.e. oscillation within a two patch length, which may not be so desirable except for when the boundary condition has a similar oscillatory nature. Secondly, even when one factor is taken as a positive value the other, obtained as the roots of the characteristic equation, can still be negative (or conjugate pairs with negative real parts). Therefore, since both ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ are taken as integral argument, negative values for one, ${\textstyle {\mu }_{1}}$ for instance, are partially covered by the values obtained form insertion of the other factor in the equation, ${\textstyle {\mu }_{2}}$ for example. Our experience shows that this form of integration is sufficient and convergent even when an oscillatory boundary condition exists and in fact it helps to reduce the number of undesirable oscillatory modes which are to be removed from the solution.

The above integrals may be taken numerically by a Gauss quadrature rule, for instance, excluding ${\textstyle \mu =0}$ from the series of integration points.

It may be noticed that ${\textstyle \vert {\mu }_{1}\vert =1}$ (or ${\textstyle \vert {\mu }_{2}\vert =1}$ ) does not correspond to decay condition but in fact periodic one. But this is of no importance because of two reasons. Firstly, if such periodic condition is not excluded from the asymptotic solution, it will not be independent of ${\textstyle {\boldsymbol {\overline {u}}}_{h}^{p}}$ initially considered for the interior of the domain and thus the formulation will inherently include both solutions, i.e. periodic ones and the correction for the boundaries. Secondly, if the boundary conditions specified at ${\textstyle x=0}$ or ${\textstyle y=0}$ inherently reflect the decay condition, i.e. vanish for large coordinates, then such a non decay solution vanishes automatically.

Supposing that the coefficient values are available we construct the finite element solution ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ for the whole domain through defining characteristic vector ${\textstyle {\boldsymbol {w}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}}$ from its smallest subset, i.e. ${\textstyle {\phi }_{{\mbox{ }}{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}}$ , using a transformation matrix like ${\textstyle {\boldsymbol {T^{'}}}_{\left({\mu }_{k},h({\mu }_{k})\right)}^{{\mbox{ }}{\mbox{ }}k}}$ containing proportionality factors

(64)

Clearly the dimension of depends upon the number of patches considered along two axes. Therefore correction to finite element solution due to presence of boundary condition is evaluated as

(65)

Now the question regarding the values of coefficients, ${\textstyle c_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}}$ , must be answered. For this we revisit the cases for boundary conditions.

1. Dirichlet boundary condition: Before imposing the boundary conditions at ${\textstyle x=0}$ or ${\textstyle y=0}$ , the nodal values along the two axes should be taken out from the vectors in null space of (Figure 4-b). This may be accomplished by defining a suitable transformation matrix as

(66)

Here ${\textstyle {\boldsymbol {T}}_{0}^{k}}$ contains the proportionality factors, ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ , and selects the nodal values of at boundaries. Here again the dimension of depends on the number of the boundary patches considered along two axes. Only few patches are practically needed for this purpose. The prescribed boundary values in this case are as (32) instead of which, for simplicity, we use ${\textstyle {\boldsymbol {\overline {u}}}_{prs}}$ . Therefore

(67)

We assume that the coefficient value associated with a value for ${\textstyle {\mu }_{k}}$ is proportional to the projection of on . Such a proportionality assumption, which is consistent with conventional mathematical transformations, is applied in most general form as

(68)

Where is a matrix to be determined so that Equation (68) holds for all values of ${\textstyle {\mu }_{k}}$ . To evaluate we substitute (68) in (67) to obtain

(69)

Since is assumed to be independent of ${\textstyle {\mu }_{k}}$ it can be taken out of the integral

{\boldsymbol {\overline {u}}}_{prs}=\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v}}_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {R}}{\mbox{ }}{\boldsymbol {\overline {u}}}_{prs}

(70)

Above relation implies that

{\boldsymbol {R}}={\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v}}_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}\right]d{\mu }_{k}\right)}^{-1}

(71)

Our experience shows that ${\textstyle {\boldsymbol {R}}}$ is well-conditioned and its inverse exists as long as enough integration points are used. The number of integration points required is dependent on the number of prescribed values, i.e. the length of ${\textstyle {\boldsymbol {\overline {u}}}_{prs}}$ . The reader will note that vectors ${\textstyle {\boldsymbol {v}}_{i}}$ are not orthogonal and thus increasing the number of integration points helps to obtain a full rank matrix.

Having found the proportionality matrix , the finite element solution corresponding to the correction due to the presence of boundary is evaluated from (65) as

{\boldsymbol {\overline {u}}}_{h0}^{p}=\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum \left\{{\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}{\boldsymbol {R}}{\boldsymbol {\overline {u}}}_{prs}\right\}{\boldsymbol {w}}_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right]{\mbox{ }}d{\mu }_{k}

(72)

or

{\boldsymbol {\overline {u}}}_{h0}^{p}=\left[\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left\{\sum {\boldsymbol {w}}_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}\right\}d{\mu }_{k}\right]{\boldsymbol {R}}{\mbox{ }}{\boldsymbol {\overline {u}}}_{prs}

(73)

The rest of the procedure is construction of general asymptotic finite element solution as mentioned before.

2. Neumann boundary condition: For this type of boundary conditions we revisit relation (37) as

(74)

where ${\textstyle {\boldsymbol {B}}_{bnd}}$ denotes a partition of ${\textstyle {\boldsymbol {B}}}$ corresponding to the set of nodes at boundary on which ${\textstyle {\boldsymbol {F}}_{prs}}$ is applied and ${\textstyle {\boldsymbol {B}}_{1}}$ is the partition pertinent to nodes the shape function of which have projection to the nodes of ${\textstyle {\boldsymbol {B}}_{bnd}}$ . Consequently ${\textstyle {\left({\boldsymbol {\overline {u}}}_{h0}^{p}\right)}_{bnd}}$ in (74) represents a partition of ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ corresponding to nodes at or near boundary. More briefly

(75)

Here ${\textstyle {\boldsymbol {K}}_{bnd}}$ is a rectangular matrix of coefficients. We note that since in (74) just elements with at least one node on boundaries are contributing, a suitable transformation matrix as below should be defined

(76)

Where ${\textstyle {\boldsymbol {T^{'}}}_{0}^{{\mbox{ }}{\mbox{ }}k}}$ contains the proportionality factors ${\textstyle {\mu }_{1}}$ and ${\textstyle {\mu }_{2}}$ and selects the nodal values of ${\boldsymbol {\phi }}_{{\mbox{ }}{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}$ at or near boundaries.

Noting that a suitable subset of ${\textstyle {\boldsymbol {w^{k}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}}$ as ${\textstyle {\boldsymbol {v^{'k}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}}$ , similar to ${\textstyle {\left({\boldsymbol {\overline {u}}}_{h0}^{p}\right)}_{bnd}}$ which is a subset of ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ , must be used then substitution of (65) in (75) leads to

(77)

or

(78)

In above ${\textstyle {\boldsymbol {K}}_{bnd}{\boldsymbol {v}}_{i({\mu }_{k})}^{\prime k}}$ represents the contribution of ${\textstyle {\boldsymbol {w^{k}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}}$ to nodal boundary resultants of tractions. Therefore, similar to the first case, we assume that the coefficient value associated with a value for ${\textstyle \mu _{k}}$ is proportional to the projection of ${\textstyle {\boldsymbol {K}}_{bnd}{\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}}$ on ${\textstyle {\boldsymbol {F}}_{prs}}$ . Hence

(79)

Where, again, ${\textstyle {\boldsymbol {R^{'}}}}$ is to be determined so that Equation (79) holds for all values of ${\textstyle {\mu }_{k}}$ . To evaluate ${\textstyle {\boldsymbol {R^{'}}}}$ we substitute (79) in (78) as

\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {K}}_{bnd}{\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}{\left({\boldsymbol {K}}_{bnd}{\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]{\boldsymbol {R^{'}}}{\mbox{ }}{\boldsymbol {F}}_{prs}d{\mu }_{k}=

(80)

and thus

{\boldsymbol {K}}_{bnd}\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}{\left({\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {K}}_{bnd}^{T}{\boldsymbol {R^{'}}}{\mbox{ }}{\boldsymbol {F}}_{prs}=

(81)

Which implies that

{\textbf {R}}'=\left[{\textbf {K}}_{bnd}\left(\sum _{k=1}^{2}\int _{0}^{1}\left[\sum {\textbf {v}}_{i(\mu _{k})}^{'k}({\textbf {v}}_{i(\mu _{k})}^{'k})^{T}\right]d\mu _{k}\right){\textbf {K}}_{bnd}^{T}\right]^{-1}

(82)

The finite element solution corresponding to the correction due to the presence of boundary is then evaluated from (65) as

{\bar {\bf {u}}}_{h0}^{p}=\sum _{k=1}^{2}\int _{\,0}^{\,1}\left(\sum {\textbf {w}}_{\,i\,(\mu _{k})}^{k}\left[\left({\textbf {K}}_{bnd}{\textbf {v}}_{i\,(\mu _{k})}^{'k}\right)^{T}{\textbf {R}}'\;{\textbf {F}}_{prs}\right]\right)d\mu _{k}

(83)

or

{\boldsymbol {\overline {u}}}_{h0}^{p}=\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {w}}_{{\mbox{ }}i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {K}}_{bnd}^{T}{\boldsymbol {R^{'}}}{\mbox{ }}{\boldsymbol {F}}_{prs}

(84)

3. General mixed boundary condition. To obtain a general relation between nodal resultants of tractions and associative nodal values, one may use (79) in (67). Therefore by dropping the subscripts for ${\textstyle {\boldsymbol {F}}}$ and we have

{\boldsymbol {\overline {u}}}=\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {K}}_{bnd}{\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]{\boldsymbol {R^{'}}}{\mbox{ }}{\boldsymbol {F}}{\mbox{ }}d{\mu }_{k}

(85)

or

(86)

in flexibility form or in stiffness form as

{\boldsymbol {F}}={\left[\sum _{k=1}^{2}\left({\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]d{\mu }_{k}\right){\mbox{ }}{\boldsymbol {K}}_{bnd}^{T}{\boldsymbol {R^{'}}}\right]}^{-1}{\boldsymbol {\overline {u}}}

(87)

Equation (87) allows us to use any mixed boundary conditions similar to a standard finite element solution. Of course the associative ${\textstyle {\boldsymbol {\overline {u}}}_{h0}^{p}}$ will be as (84) while dropping the subscript of ${\textstyle {\boldsymbol {F}}}$ .

It may be noticed that we could also use (73) in (74) for relating ${\textstyle {\boldsymbol {F}}}$ to ${\textstyle {\boldsymbol {\overline {u}}}}$ in one set of equations. In that case ${\textstyle {\boldsymbol {w}}}$ should be replaced by ${\textstyle {\boldsymbol {v^{'}}}}$ and therefore;

{\boldsymbol {K}}_{bnd}\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[{\sum {\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {R}}{\mbox{ }}{\boldsymbol {\overline {u}}}=

(88)

It follows that theoretically the following relation must exist between and

{\boldsymbol {K}}_{bnd}\left(\sum _{k=1}^{2}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[{\sum {\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\left({\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {R}}=

{\left[\sum _{k=1}^{2}\left({\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\left[\sum {\boldsymbol {v}}_{i{\mbox{ }}{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{k}{\left({\boldsymbol {v^{'}}}_{i{\mbox{ }}{\mbox{ }}({\mu }_{k})}^{{\mbox{ }}k}\right)}^{T}\right]d{\mu }_{k}\right){\boldsymbol {K}}_{bnd}^{T}{\boldsymbol {R^{'}}}\right]}^{-1}

(89)

Of course from a numerical point of view there may be some discrepancy between the two relations.

It is noteworthy to mention that the boundary conditions specified should be consistent with assumption made for the decay property, i.e. ${\textstyle {\boldsymbol {\overline {u}}}_{prs}\rightarrow 0}$ or ${\textstyle {\boldsymbol {F}}_{prs}\rightarrow 0}$ when ${\textstyle {\frac {d_{x}}{L_{x}}}\rightarrow \infty }$ and ${\textstyle {\frac {d_{y}}{L_{y}}}\rightarrow \infty }$ which of course is the case for our studies.

If asymptotic finite element solution near a corner of a domain is sought, the decay condition ${\textstyle {\frac {d_{z}}{L_{z}}}\rightarrow \infty }$ is also added to the above conditions. To satisfy such an additional decay condition a corresponding proportionality property, with ${\textstyle {\mu }_{3}}$ as a parameter, between nodal values along ${\textstyle z}$ direction must be considered. The procedure leads to finding a characteristic equation with three parameters

(90)

form which one parameter is determined by assuming values for the other two. We define


${\mu }_{3}=h({\mu }_{1},{\mu }_{2})$	$h=\left\{g{\mbox{ }}{\mbox{ }}{\mbox{ }}\vert {\mbox{ }}{\mbox{ }}f({\mu }_{1},{\mu }_{2},g)=\right.$ $\left.0,{\mbox{ }}{\mbox{ }}{\mbox{ }}\vert g\vert \leq 1,{\mbox{ }}{\mbox{ }}\vert {\mu }_{1}\vert \leq 1,{\mbox{ }}{\mbox{ }}\vert {\mu }_{2}\vert \leq 1\right\}$

| style="width: 5px;text-align: right;white-space: nowrap;" | (91) |}

Similar expressions may be written by replacing the parameters in turn. Having found the null space of ${\textstyle {\boldsymbol {Q^{'}}}}$ for each set, the independent nodal values of the patch is evaluated as

(92)

or

(93)

We have omitted the superscripts for convenience. Finally the formulation leads to

{\boldsymbol {\overline {u}}}_{h0}^{p}=\left[\sum _{k=1}^{2}\sum _{n=k+1}^{3}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}{\int }_{{\mbox{ }}0}^{{\mbox{ }}1}\sum {\boldsymbol {w}}_{i{\mbox{ }}({\mu }_{k},{\mu }_{n})}{\boldsymbol {v}}_{i{\mbox{ }}({\mu }_{k},{\mu }_{n})}^{T}d{\mu }_{n}d{\mu }_{k}\right]{\boldsymbol {R}}{\mbox{ }}{\boldsymbol {\overline {u}}}_{prs}

(94)

Where

{\boldsymbol {R}}={\left[\sum _{k=1}^{2}\sum _{n=k+1}^{3}{\int }_{0}^{1}{\int }_{0}^{1}\sum {\boldsymbol {v}}_{i{\mbox{ }}({\mu }_{k},{\mu }_{n})}{\boldsymbol {v}}_{i{\mbox{ }}({\mu }_{k},{\mu }_{n})}^{T}d{\mu }_{n}d{\mu }_{k}\right]}^{-1}

(95)

Other relations may be found by analogy.

It is observed that a finite element solution is achievable for Equation (49) with boundary conditions specified in (50) and (51). The accuracy obtained depends upon the number of boundary patches nominated for satisfaction of boundary conditions at ${\textstyle x=0}$ , ${\textstyle y=0}$ and ${\textstyle z=0}$ . In the limit, when the number of boundary patches tends to infinity, the solution must approach to a unique solution since the integral equation (49) and all boundary conditions are satisfied in a finite element sense. Our experience shows that few boundary patches are needed to obtain a reasonable converged solution. In a following sample problem we examine the convergence of the solution for a two dimensional case.

4. NUMERICAL PRESENTATION

A sample problem

Apart from assessment of error estimators, the method described above may be used to find a solution to a homogeneous differential equation like following Poisson’s equation

and

with exact solution, for instance, as

${\textstyle u_{ex}=e^{-(x+y)}cos(x-y)}$

It may be seen that the decay condition exists in the exact solution when ${\textstyle x\rightarrow \infty }$ or ${\textstyle y\rightarrow \infty }$ . We try to give a finite element solution to such a problem using the proposed approach.


(a)	(b)	(c)

Figure 5. Patch used for sample problem; (a) Union Jack pattern, (b) Four square elements as the basic patch, (c) A square element as the basic patch.

Solution with triangular elements

A mesh of infinite dimensions is considered with a repeatable pattern in which the smallest unit is of form as the patch shown in Figure 5-a. After assembling the stiffness matrix of a macro patch of ${\textstyle 2\times 2}$ from the smallest unit, condensing the matrix for independent nodal values and also imposing the proportionality conditions similar to relations (52), the following expression for ${\textstyle {\boldsymbol {Q^{'}}}}$ is obtained

${\boldsymbol {Q^{'}}}\left({\mu }_{1},{\mu }_{2}\right)=$ $\left[{\begin{array}{ccc}4&-1-{\frac {1}{{\mu }_{1}}}&-1-{\frac {1}{{\mu }_{2}}}\\-1-{\mu }_{1}&{\frac {7}{2}}-{\frac {1}{4{\mu }_{2}}}-{\frac {{\mu }_{2}}{4}}&-{\frac {1}{4}}-{\frac {{\mu }_{1}}{4}}-{\frac {1}{4{\mu }_{2}}}-{\frac {{\mu }_{1}}{4{\mu }_{2}}}\\-1-{\mu }_{2}&-{\frac {1}{4}}-{\frac {1}{4{\mu }_{1}}}-{\frac {{\mu }_{2}}{4}}-{\frac {{\mu }_{2}}{4{\mu }_{1}}}&{\frac {7}{2}}-{\frac {1}{4{\mu }_{1}}}-{\frac {{\mu }_{1}}{4}}\end{array}}\right]$

The determinant of ${\textstyle {\boldsymbol {Q^{'}}}}$ is set to be zero and thus the characteristic equation is obtained as

$f\left({\mu }_{1},{\mu }_{2}\right)=33+{\frac {1}{4}}\left({\mu }_{1}^{2}+\right.$ $\left.{\mu }_{2}^{2}+{\mu }_{1}^{-2}+{\mu }_{2}^{-2}\right)-$ $8\left({\mu }_{1}+{\mu }_{2}+{\mu }_{1}^{-1}+{\mu }_{2}^{-1}\right)-$ ${\frac {1}{2}}\left({\mu }_{1}{\mu }_{2}+{\mu }_{1}^{-1}{\mu }_{2}+\right.$ $\left.{\mu }_{1}{\mu }_{2}^{-1}+{\mu }_{1}^{-1}{\mu }_{2}^{-1}\right)=$ $0$

We evaluate ${\textstyle {\mu }_{2}}$ in terms of ${\textstyle {\mu }_{1}}$

{\mu }_{2}=g_{1}({\mu }_{1})={\frac {1}{2}}\left(16+{\mu }_{1}^{-1}-8{\mu }_{1}^{-{\frac {1}{2}}}-8{\mu }_{1}^{\frac {1}{2}}+{\mu }_{1}-{\mu }_{1}^{-1}\left(-1+5{\mu }_{1}^{\frac {1}{2}}-5{\mu }_{1}+{\mu }_{1}^{\frac {3}{2}}\right){\left(1-6{\mu }_{1}^{\frac {1}{2}}+{\mu }_{1}\right)}^{\frac {1}{2}}\right)

${\mu }_{2}=g_{2}\left({\mu }_{1}\right)={\frac {1}{2}}\left(16+{\mu }_{1}^{-1}+8{\mu }_{1}^{-{\frac {1}{2}}}+8{\mu }_{1}^{\frac {1}{2}}+{\mu }_{1}-{\mu }_{1}^{-1}\left(1+5{\mu }_{1}^{\frac {1}{2}}+5{\mu }_{1}+{\mu }_{1}^{\frac {3}{2}}\right){\left(1+6{\mu }_{1}^{\frac {1}{2}}+{\mu }_{1}\right)}^{\frac {1}{2}}\right)$

${\mu }_{2}=g_{3}\left({\mu }_{1}\right)=g_{1}{\left({\mu }_{1}\right)}^{-1},{\mbox{ }}{\mu }_{2}=g_{4}\left({\mu }_{2}\right)=g_{2}{\left({\mu }_{1}\right)}^{-1}$

Similar expressions may be derived for ${\textstyle {\mu }_{1}}$ in terms of ${\textstyle {\mu }_{2}}$ . Generally, such explicit relations are not necessary and what needed is just choosing a value for ${\textstyle {\mu }_{1}}$ , as an integration point in the feasible domain, and solving an algebraic equation for its roots. Nevertheless, to give some insight to the problem the explicit relations are given here. Now we define

,

when ${\textstyle \vert {\mu }_{1}\vert \leq 1}$ . Similar functions are defined, with analogy, when ${\textstyle {\mu }_{1}}$ is evaluated in terms of ${\textstyle {\mu }_{2}}$ .

We should note that generally all quantities may take complex values, however, as pointed out before we take real values for one parameter and calculate the other one which may be in the form of a complex pair.

To give a finite element solution, we take exact fluxes normal to boundaries at ${\textstyle x=0}$ and ${\textstyle y=0}$ as tractions for the infinite domain to be solved. We note that although relations (87) or (88) may be employed to first find the unknown nodal boundary values and then extend the result to the interior parts through Equation (73) , one may directly use Equation (84) instead. Figure 6 depicts the variation of the finite element solution for ${\textstyle u}$ obtained from the proposed method using 10 patches along each axis, ${\textstyle x=0}$ and ${\textstyle y=0}$ . In order to show that the solution is close to an approximate finite element solution, in Figure 7 we give the results of a finite element solution on a domain of ${\textstyle 10\times 10}$ patches using exact fluxes for the other two boundaries. Although we do not expect that the two sets of results be identical, due to the fact that the tractions used for the latter finite element solution do not represent the tractions of remaining domain in FEM sense, the results are very close (the reader may also note that since all the boundary conditions for the latter solution are in Neumann form one node of the domain should be restrained, the most remote one in this case, thus this constitutes another source for the difference between the two sets of results). It can clearly be seen that the method proposed is able to give a finite element solution for a domain with unbounded dimension exhibiting a decay property.


(a)	(b)

Figure 6. Finite solution obtained from the proposed method for the sample problem using ${\textstyle 10\times 10}$ patches with Union Jack pattern in Fig. 5-a and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.


(a)	(b)

Figure 7. Standard finite solution for the sample problem using ${\textstyle 10\times 10}$ patches with Union Jack pattern in Fig. 5-a and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ , and other two boundaries; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.

The accuracy of the solution may be increased or decreased by use of different number of patches along the two axes. Figures 8 and 9 show the variation of ${\textstyle u}$ over the domain using different number of patches along each direction.


(a)	(b)

Figure 8. Finite solution obtained from the proposed method for the sample problem using ${\textstyle 8\times 8}$ patches with Union Jack pattern in Fig. 5-a and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ . The results are shown on ${\textstyle 10\times 10}$ patches; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.


(a)	(b)

Figure 9. Finite solution obtained from the proposed method for the sample problem using ${\textstyle 6\times 6}$ patches with Union Jack pattern in Fig. 5-a and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ . The results are shown on ${\textstyle 10\times 10}$ patches; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.

Figure 10 also show the convergence of the results in terms of the number of integration points and patches used. In these figures the nodal values obtained for nodes 1 and 5 of the corner patch in Figure 4-b are plotted against the number of Gauss points and patches used. At the first glimpse no uniform convergence is seen, however, more careful look into the results reveals that the two values for nodes 1 and 5 follow each other with constant difference. The reason for such behavior lies in the rigid body motion of the whole sub-domain. We notice that the boundary conditions used are of Neumann type and no explicit essential condition is used except for the decay condition the strength of which varies between different solutions. The gradients and relative values show excellent convergence for moderate number of Gauss points and patches as shown in Figure 10 for ${\textstyle u_{5}-u_{1}}$ and ${\textstyle {\sigma }_{x}}$ , the gradient at center of element 1-2-9 in Figure 5-a.


(a)	(b)

Figure 10. Convergence of the solution with Neumann boundary condition- ${\textstyle u_{1}}$ and ${\textstyle u_{5}}$ are nodal values of node one and five in Union Jack Pattern of Fig 5-a for the corner patch of the domain, and ${\textstyle \sigma }$ is the gradient at center of triangle 1-2-9; (a) Convergence with respect to the number of Gauss points in a ${\textstyle 10\times 10}$ macro patch, (b) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.

We now examine the performance of the method when Dirichlet type of boundary conditions is of concern. The exact values are taken as prescribed boundary conditions along ${\textstyle x=0}$ and ${\textstyle y=0}$ . Equation (73) is directly employed to obtain the solution. To give a finite element solution for comparison, here again a square domain is solved with the prescribed values, at ${\textstyle x=0}$ and ${\textstyle y=0}$ , and Neumann boundary conditions at the other two edges. The results of the two solutions are presented in Figures 11 and 12. It can again be seen that the proposed method is able to produce the finite element solution. Figure 13 depicts the convergence of the proposed solution, in terms of the number of integration points and patches used. Unlike the case for Neumann boundary conditions, in this case uniform and fast convergence is seen in the results.


(a)	(b)

Figure 11. Finite solution obtained from the proposed method for the sample problem using ${\textstyle 10\times 10}$ patches with Union Jack pattern in Fig. 5-a and Dirichlet boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.


(a)	(b)

Figure 12. Standard finite solution for the sample problem using ${\textstyle 10\times 10}$ patches with Union Jack pattern in Fig. 5-a and Dirichlet boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ , and Neumann condition at the other two boundaries; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.


(a)	(b)

Figure 13. Convergence of the solution with Dirichlet boundary condition- ${\textstyle u_{1}}$ and ${\textstyle u_{5}}$ are nodal values of node one and five in Union Jack Pattern of Fig 5-a for the corner patch of the domain, and ${\textstyle \sigma }$ is the gradient at center of triangle 1-2-9; (a) Convergence with respect to the number of Gauss points in a ${\textstyle 10\times 10}$ macro patch, (b) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.

Solution with four node elements

A mesh consisting of repeated patterns of a four square element as a macro-patch shown in Figure 5-b is considered for solution. Applying the proportionality relations, as (52), the following expression for the characteristic equation is obtained

$f\left({\mu }_{1},{\mu }_{2}\right)={\frac {799}{54}}-{\frac {259}{108}}F_{1}+$ ${\frac {25}{54}}F_{2}-{\frac {1}{108}}F_{3}-{\frac {49}{108}}F_{1}^{2}+$ ${\frac {1}{108}}F_{1}F_{2}-{\frac {1}{648}}F_{1}F_{3}+{\frac {1}{864}}F_{2}^{2}+$ ${\frac {1}{2592}}F_{1}^{4}$

Where

$F_{1}={\mu }_{1}+{\mu }_{1}^{-1}+{\mu }_{2}+{\mu }_{2}^{-1}$ , $F_{2}={\mu }_{1}^{2}+{\mu }_{1}^{-2}+{\mu }_{2}^{2}+{\mu }_{2}^{-2}$ , $F_{3}={\mu }_{1}^{3}+{\mu }_{1}^{-3}+{\mu }_{2}^{3}+{\mu }_{2}^{-3}$

From above characteristic equation ${\textstyle {\mu }_{2}}$ is to be calculated in terms of values given for ${\textstyle {\mu }_{1}}$ and vice versa. Two forms of problem, i.e. with Neumann and Dirichlet conditions, are solved using the proposed method. Some results obtained from using different numbers of patches along the axes are given in Figures 14 and 15. Once again, in order to compare the results with a standard finite element solution, a square domain is solved with the conditions mentioned previously, i.e. tractions from exact solution at all edges (see Figure 16). Convergence of the solution for one node of the domain is illustrated in Figure 17 in terms of number of integration points and patches used. Similar discussion as given for Figure 10 is valid here and as is seen the relative quantities and the gradients exhibit excellent convergence.


(a)	(b)

Figure 14. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of four squares in Fig. 5-b and Dirichlet boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Solution for ${\textstyle 8\times 8}$ patches shown on a domain of ${\textstyle 10\times 10}$ patches, (b) Solution for ${\textstyle 10\times 10}$ patches.


(a)	(b)

Figure 15. Contour plots the finite solutions obtained from the proposed method for the sample problem using pattern of four squares in Fig. 5-b and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Solution for ${\textstyle 8\times 8}$ patches shown on a domain of ${\textstyle 10\times 10}$ patches, (b) Solution for ${\textstyle 10\times 10}$ patches.


(a)	(b)

Figure 16. Contour plots for standard finite solutions for the sample problem using ${\textstyle 10\times 10}$ pattern of four squares in Fig. 5-b and boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Neumann condition, (b) Dirichlet condition.


(a)	(b)

(c)	(d)

Figure 17. Convergence of the solution with Neumann and Dirichlet boundary condition- ${\textstyle u_{1}}$ and ${\textstyle u_{5}}$ are nodal values of node one and five in Fig 5-b for the corner patch of the domain, and ${\textstyle \sigma }$ is the gradient at center of element 1-2-9-8; (a)/(c) Convergence of solution with respect to the number of Gauss points in a ${\textstyle 10\times 10}$ macro patch, (b)/(d) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.

It is interesting to note that for solution using square elements, the smallest repeatable pattern may be considered as one element. Therefore, it is worthwhile to compare the results from such a solution with those obtained above.

If the proportionality conditions are written for one element, just one relation is obtained for ${\textstyle {\boldsymbol {Q^{'}}}}$ which is identical to the characteristic equation

$f\left({\mu }_{1},{\mu }_{2}\right)={\frac {1}{3}}\left[8-\right.$ $\left.\left({\mu }_{1}+{\mu }_{1}^{-1}+{\mu }_{2}+{\mu }_{2}^{-1}+\right.\right.$ $\left.\left.{\mu }_{1}{\mu }_{2}+{\mu }_{1}{\mu }_{2}^{-1}+\right.\right.$ $\left.\left.{\mu }_{1}^{-1}{\mu }_{2}+{\mu }_{1}^{-1}{\mu }_{2}^{-1}\right)\right]$

Figures 18 to 20 show the results obtained. Here ${\textstyle 20\times 20}$ patches are used in order that the results be comparable with those of Figures 14,15 and 17. Study on the convergence of the solution, while comparing to the solution with ${\textstyle 2\times 2}$ elements, demonstrates a slower rate of convergence in terms of the number of Gauss points and a faster rate in terms of the number of patches used (Figure 20).


(a)	(b)

Figure 18. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of one square element in Fig. 5-c and Dirichlet boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Solution for ${\textstyle 16\times 16}$ patches shown on a domain of ${\textstyle 20\times 20}$ patches, (b) Solution for ${\textstyle 20\times 20}$ patches.


(a)	(b)

Figure 19. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of one square element in Fig. 5-c and Neumann boundary condition at ${\textstyle x=0}$ and ${\textstyle y=0}$ ; (a) Solution for ${\textstyle 16\times 16}$ patches shown on a domain of ${\textstyle 20\times 20}$ patches, (b) Solution for ${\textstyle 20\times 20}$ patches.


(a)	(b)

(c)	(d)

Figure 20. Convergence of the solution with Neumann and Dirichlet boundary condition- ${\textstyle u_{1}}$ and ${\textstyle u_{5}}$ are nodal values of node one and five in Fig 5-c for the corner patch of the domain, and ${\textstyle \sigma }$ is the gradient at center of element 1-2-3-4; (a)/(c) Convergence of solution with respect to the number of Gauss points in a ${\textstyle 20\times 20}$ macro patch, (b)/(d) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.

5. CONCLUSIONS

The formulation introduced by Babuška and co workers [13], for assessment of error estimator near flat boundaries, has been extended to kinked boundaries. The new formulation given is suitable for study of either asymptotic behavior of a finite element solution near boundaries or general FEM solution for unbounded domains.

Asymptotic finite element solution consists of two main parts; one obtained from patches located at interior of the domain, and another obtained from a boundary layer solution irrespective of other boundaries at far ends of the domain. The latter part of the solution decays towards the interior of the domain and, in this respect, is similar to a finite element solution of an unbounded domain.

The proposed method is based on assuming proportionality relations between the finite element nodal values along two or three directions, when the mesh is of a repeatable pattern. Writing the equilibrium equations over a unit patch leads to a characteristic equation containing all proportionality factors from which one may be expressed in terms of the others. The boundary layer solution is obtained by integral of extended vectors arranged in spectral form from null space vectors of the equilibrium equations. The decay condition is then applied by defining a feasible domain for proportionality factors.

Performance of the method has been examined for solution of unbounded domains in order to ensure its applicability to producing asymptotic finite element solution in an assessment procedure. A sample problem has been solved using two forms of patterns to demonstrate the capability of the method in reproducing the boundary layer solution. The convergence of the procedure has been studied with respect to number of patches as well as the number of integration points used. The results show fast convergence for both cases.

In the next part of the paper we shall give the step by step procedure for robustness test followed by numerous results obtained from application of the test to some recovery based error estimators.

APPENDIX

In this appendix we explain the reason for existence of proportionality property used in (52-a) to (52-c).

Integral Equation (27) may be viewed as a finite element solution of a corresponding differential equation with constant coefficient as the following one

on

(A-1)

For scalar problems such as heat conduction ones, for instance, the differential equation becomes

(A-2)

forgetting the subscript and superscripts for simplicity. Since the equation is of constant coefficient type, the solution may be expressed in an exponential form as ${\textstyle u=ce^{\alpha x+\beta y}}$ . This leads to a characteristic equation as ${\textstyle {\alpha }^{2}+{\beta }^{2}=0}$ after substitution in the equation. The exponential form of solution of (A-2) may be written as (see also [17])

(A-3)

or

(A-4)

where integrations are taken on feasible domain of ${\textstyle \alpha }$ or ${\textstyle \beta }$ , e.g. conditions for decay property. Considering single solution from (A-3) for example, it may be seen that the following relations exist

Similarly

or more generally

$u_{\left(x+nL_{x},y+mL_{y}\right)}=c_{\left(\alpha \right)}e^{\alpha (x+nL_{x})+\beta (y+mL_{y})}=$ ${\left({\mu }_{01}\right)}^{n}{\left({\mu }_{02}\right)}^{m}u_{\left(x,y\right)}$

Which is consistent with assumption made in (52-a) to (52-c). We note that here the differential equation is satisfied exactly whereas in relations (52) the aim is to satisfy the integral form of the equation. Of course the proportionality factors in the two approaches are not similar since one is an approximate finite element solution of another.

An important observation may be made by looking in to the characteristic equation ${\textstyle {\alpha }^{2}+{\beta }^{2}=0}$ . It can be seen that for a given value of ${\textstyle \alpha }$ there will be two values for ${\textstyle \beta }$ one of which corresponds to decay condition, i.e. ${\textstyle e^{-i\alpha y}}$ or ${\textstyle e^{i\alpha y}}$ (note that ${\textstyle \alpha }$ can be a complex value). Hence for a given value of ${\textstyle {\mu }_{01}}$ , two values as ${\textstyle {\mu }_{02}}$ and ${\textstyle {\mu }_{02}^{-1}}$ are both feasible but the one with absolute value less than unity corresponds to decay condition. Similar relations may be written for a general three dimensional problem.

References

[1] Oden, J.T., Belytschko, T., Babuška, I. and Hughes, T.J.R., “Research directions in computational mechanics” Comput. Methods Appl. Mech. Eng., Vol. 192, pp. 913-922, 2003.

[2] Zienkiewicz, O.C., “Achievements and some unsolved problems of finite element method” , Int. J. Numer. Meth. Engrg., Vol. 47, pp. 9-28, 2000.

[3] Hinton, E. and Campbell, J.S., “Local and global smoothing of discontinuous finite element functions using a least square method”, Int. J. Numer. Meth. Engrg., Vol. 8, pp. 461-480, 1976.

[4] Zienkiewicz, O.C. and Zhu, Z., “A simple error estimator and adaptivity procedure for practical engineering analysis”, Int. J. Numer. Meth. Engrg., Vol. 24, pp. 337-357, 1987.

[5] Tessler, A., Riggs, H.R. and Dambach, M., “A novel four-node quadrilateral smoothing element for stress enhancement and error estimation”, Int. J. Numer. Meth. Engrg., Vol. 44, pp. 1527-1543, 1999.

[6] Hiller, J.F. and Bathe, K.J., “On higher-order-accuracy points in isotropic finite element analysis and application to error assessment”, Comput. & Struc., Vol. 79, pp. 1275-1285, 2001.

[7] Zienkiewicz, O.C. and Zhu, Z., “The superconvergent patch recovery and a posteriori error estimates. Part II: Error estimates and adaptivity”, Int. J. Numer. Meth. Engrg., Vol. 33, pp. 1365-1382, 1992.

[8] Babuška, I., Strouboulis, T. and Upadhyay, C.S., “ A model study of the quality of a posteriori error estimaors for linear elliptic problems. Error estimation interior of patchwise uniform grids of triangles”, Comput. Methods Appl. Mech. Eng., Vol. 114, pp. 307-378, 1994.

[9] Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K., “Validation of a posteriori error estimators by numerical approach”, Int. J. Numer. Meth. Engrg., Vol. 37, pp. 1073-1123, 1994.

[10] Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K. “An objective criterion for assessing the reliability of a posteriori error estimators in finite element computations”, U.S.A.C.M. Bulletin, No.7 , pp. 4-16, 1994.

[11] Zhang, Z. and Zhang, S., “ Derivative superconvergence of rectangular finite elements for the Reissner-Mindlin plate”, Comput. Methods Appl. Mech. Eng., Vol. 134, pp. 1-16, 1996.

[12] Zhang, Z., “ Superconvergence in projected –shear plate-bending finite element methods”, Numer. Meth. Partial Diff. Eqn., Vol. 14, pp. 367-386, 1998.

[13] Babuška, I., Strouboulis, T. and Upadhyay, C.S., “A model study of the quality of a posteriori error estimators for finite element solutions of linear elliptic problems, with particular reference to the behavior near the boundary”, Int. J. Numer. Meth. Engrg., Vol. 40, pp. 2521-2577, 1997.

[14] Ainsworth, M. and Oden, J.T., A posteriori error estimation in finite element analysis, John Wiley, 2000.

[15] Zienkiewicz, O. C. and Taylor, R. L., “The finite element method”, 5^th edition, Butterworth-Heineman, 2000.

[16] Boroomand, B. and Zienkiewicz, O.C., “An improved REP recovery and efectivity robustness test”, Int. J. Numer. Meth. Engrg., Vol. 33, pp. 3247-3277, 1997.

[17] Hildebrand, F., Advanced calculus for applications, Second edition, Prentice-Hall, 1976.

Latest revision as of 12:50, 15 June 2018

Abstract

1. INTRODUCTION

2. PRELIMINARIES

3. ROBUSTNESS PATCH TEST

4. NUMERICAL PRESENTATION

5. CONCLUSIONS

APPENDIX

References

Document information

Document Score

Share this document

Keywords

claim authorship

@@ Line 1: / Line 1: @@
+==Abstract==
-==On application of robustness test for error estimators to patches near kinked boundaries in three dimensional problems: Part II,  Test results for error estimators using SPR and REP==
+In this part of paper we shall extend the formulation proposed by Babuška and co-workers for robustness patch test, in quality assessment of error estimators, to more general cases of patch locations especially in three dimensional problems.  This is performed first by finding an asymptotic finite element solution at interior parts of a problem with assumed smooth exact solution and then adding a correction part to obtain the solution near a kinked boundary irrespective of other boundary conditions at far ends of the domain.  It has been shown that the solution corresponding to the correction part may be obtained in an spectral form by assuming a suitable proportionality relation between the nodal values of a mesh with repeatable pattern of macro-patches.
-'''B. Boroomand''', '''F. Mossaiby'''
+Having found the asymptotic finite element solution, the performance of error estimators may be examined. Although in this paper we focus on the asymptotic behavior of error estimators, the method described in this part may be used to obtain finite element solution for two/three dimensional unbounded heat/elasticity problems with homogeneous differential equations.  Some numerical results are presented to show the validity and performance of the proposed method.
-Civil Engineering Department, Isfahan University of Technology, Iran
+'''Keywords:''' ''Recovery, Error Estimate, Superconvergent, Robustness test, three dimensional, Unbounded domains''
-==Abstract==
+=1. INTRODUCTION=
-In this part of the paper we shall use the formulation given in the first part to assess the quality of recovery based error estimators using two recovery methods, i.e. Superconvergent Patch Recovery (SPR) and Recovery by Equilibrium in Patches (REP).  The recovery methods have been shown to be asymptotically robust and superconvergent when applied to two dimensional problems. In this study we shall examine the behavior of the recovery methods on several three dimensional mesh patterns for patches located either inside or at boundaries. This is performed by first finding an asymptotic finite element solution, irrespective of boundary conditions at far ends of the domain, and then applying the recovery methods.  The test procedure near kinked boundaries is explained in a step by step manner.  The results are given in a series of tables and figures for various cases of three dimensional mesh patterns.  It has been experienced that the full superconvergent property is generally lost due to presence of boundary layer solution and the definition of the recoveries near boundaries though the results of the robustness test is still within an acceptable range.
+Efforts by pioneering engineers in late twentieth century led to understanding the necessity of error evaluation and use of adaptive procedures in finite element solutions.  Nowadays, most of commercial finite element programs are equipped with tools for error estimation and mesh generation.  Nevertheless, the need for further investigation is still felt by scientists [1].
-'''Keywords:''' ''Recovery, error estimate, superconvergent, robustness test, three dimensional''
+The knowledge of error estimation in two dimensional problems has now reached to a level of maturity [2].  However, advances in engineering and technology have led to increase of demands for three dimensional analyses and thus it is of great interest to know whether the error analyses still work well in three dimensions.  Although many of the concepts introduced in two dimensional problems may  readily be applied to three dimensions, there are many especial cases that may not comprehensively be covered by the counterparts studied in two dimensions.  Quality assessment of error estimators may be considered as one of the cases.
+There are two basic approaches employed in studies for assessment of error estimators.  The most classic approach, used almost since the first thoughts of FEM, is benchmarking. This category of approach was employed in many early studies on error estimation, see [3,4] for example, and it is still being used as a tool for giving some evidences of good performance, e.g. [5,6].  Although benchmark problems in two and three dimensional spaces have much in common, dealing with the three dimensional ones may need some more elaboration.  For example behavior of error estimators near an edge in 3D problem, where two intersecting boundary exist, may not easily be understood from the results obtained at edges of a 2D benchmark.  In fact an edge in 3D problem has much in common with a corner in a 2D problem, yet not quite similar.
-==1. Introduction==
+The second approach is employed for studies on quality of approximate error in asymptotic manner, i.e. when the element sizes tend to zero, and has roots in basic assumptions of finite element method.   Although this type of assessment approach was inherently used by some researches to give proofs for theorems and lemmas regarding the convergence of the proposed methods, e.g. see [7], for the first time a systematic computer based type of the approach was introduced by Babuška ''et al''  in a series of papers [8-10].  Other researchers followed some similar approaches to explore superconvergence in plate problems [11, 12].
-In the first part of this paper [1], we introduced a method for capturing boundary layer solution near kinked boundaries.  As we described, the boundary layer solution is a part of asymptotic finite element solution which is used in the robustness test.  In this part of the paper we shall employ the proposed method to carry out tests on some error estimates.
+The method was basically introduced for assessment of error estimators applied at interior parts of a fictitious unbounded domain and has sound mathematical bases.  The basic feature of the method is using meshes with repeatable patterns of patches in order to reduce the size of the space needed for the finite element solution.  The test procedure is performed on a solution found for a single patch and thus, in this respect, resembles to the patch tests usually used in element technology.  Since the reduction of the size of the domain, i.e. from an infinite size to a patch with finite size, is permitted under some assumptions for asymptotic behavior of the exact solution, the test results are interpreted in an asymptotic point of view.  The formulation given for two dimensional problems is general and may be applied to three dimensional problems with no much effort.
-Error estimators are usually classified into two categories, i.e. residual based error estimators and recovery based error estimator.  The first thoughts for error evaluation were in the form of the former category and the well known paper by Babuška and Rheinboldt [2] appears to be the pioneering one.  About a decade later the latter approach was introduced by Zienkiewicz and Zhu [3].  A projection method was suggested for evaluation of continuous gradient fields.  Soon it became clear that the projection method works for certain class of elements.  The reason was detected as the lack of adequate convergence of the solution.  Few years later, the authors introduced another projection method called Superconvergent Patch Recovery (SPR) and showed its performance through some benchmark problems [4].  It has also been proven that if a projection method exhibits a slightly more convergence than the finite element solution the adaptive procedure leads the solution towards the exact solution [5]. The recovery method uses information on some sampling points at element level.   The error estimator using SPR method is sometimes called “ZZ” estimator.  Many researchers soon began to use the method and many reports were published concerning the effectiveness of the method most of which admitting its superiority. Some amendments were also introduced for further increase of the convergence [6, 7] or improvement of the results near boundaries [8].
+The method was further extended to patch locations near boundaries in two dimensional problems [13].  The new formulation includes a boundary layer solution, in an asymptotic manner, for a fictitious unbounded domain.  However, the modified version of the formulation can only be applied to a domain with single flat boundary.  This means that the method may not be applied for patches near corners in two dimensional cases. Of course, this is of not much importance since the frequency of dealing with such situations is very small compared to other cases studied in two dimensional cases.
-Despite the efforts made for improvement of SPR, the original form is still the most effective and economical version.  This has been confirmed in studies performed by Bubuška ''et al'' reported in a series of papers [9-11].  One of the outcomes of the studies was a routine for testing the error estimates in an asymptotic sense which we explained in the first part of the paper [1].  Some of the error estimators studied by Bubuška and co workers are of residual based type.  It has been shown that among all error estimators suggested to the date, for problems with smooth exact solution, the recovery based one using SPR is superior in performance.
+For three dimensional problems, however, the modified version can only be used for patch locations near a flat boundary.  This means that the cases in which the patch is located near two intersecting flat boundaries, i.e. an edge, cannot be studied.   It is obvious that frequency of such cases is not comparable with that of corners in two dimensional problems.  However, if the problem is solved for edges of a three dimensional body, the formulation may also be applied for corners of a two dimensional problem.
-The success and popularity of SPR is due to two important features, i.e. simplicity and accuracy.   However, in terms of accuracy there is an argument when singularity in solution is of concern.  In other studies by  Bubuška ''et al'' it is shown that another source of error may encounter to a finite element solution called as “pollution error” [12-14].  This sort of error cannot be captured by error estimators that use local information such as those use recovery methods. This is much pronounced when the element sizes are reasonably large especially around the singular points. This is considered as a mathematical deficiency for recovery based error estimators such as “ZZ”.   However, from engineering point of view there is an argument for such deficiency which still keeps the recovery based error estimator on top.  It has been experienced that when using “ZZ” estimator in problems with singularities, the locations of the singular points are spotted very fast.  The reason for such effect lies in the fact that the error estimated near such points is usually large and thus in a real adaptive procedure the subsequent meshes become extensively fine around the singular points.  This leads to reduction of the pollution error in other parts of the domain and therefore the accuracy of the recovery results increases.  Nevertheless, according to the studies in [12-14] the criterion for termination of the adaptive procedure must carefully be revised.
+In this paper we shall extend the formulation to the situations that the patch is near an edge/corner of the domain.  We shall employ a spectral formulation in which a proportionality property plays a key role.  It will be shown that the proportionality property has roots in the nature of the governing differential equation.   As will be shown later in this paper, the governing equation for boundary layer solution in integral form is of a homogeneous type. The formulation given in this paper is not only useful for finding asymptotic boundary layer solution but also applicable to many other unbounded problems with governing equations of similar type.  This is an important feature of the method. Having found the formulation for the boundary layer solution near kinked boundaries, the results will be verified by a sample problem with exact solution on an unbounded domain.  The test results and the discussion for quality of the error estimates are given in the next part of the paper.
-Based on above discussion, it may be concluded that assumption made in [1], initially used in [9-11], will be valid for wide range of problems when the study is performed in an asymptotic sense.  In present study we focus on performance of SPR basically in three dimensional problems.  Meanwhile we shall study the behavior of another recovery method called Recovery by Equilibrium in Patches introduced in [15].  The recovery uses a weighted form of equilibrium equations; similar to FEM formulation, on patches of elements and exhibits similar behavior to SPR in two dimensional benchmarks.  The method does not need any information at special points.  The improved version of the method was introduced in [16] where the results of robustness test were given.  Since there is no need for especial information, REP may be applied to nonlinear problems such as those involving plastic deformation [17, 18]. A hybrid recovery from combination of SPR and REP can be seen in the literature [19] which may be considered as an improved version of SPR since it still needs information at sampling points and thus does not reflect full advantages of REP.
+The layout of the paper is as follows; In Section 2 we shall prepare some preliminary formulations and definitions to be used in subsequent sections.  In Section 3 some basic assumptions initially made in [8-10] and basic relations for the test at interior parts will be revisited.  The relations will be frequently referred to especially when we describe the step by step test procedure in its generalized form.  Having recalled the basic relations, the formulations used for flat boundaries, introduced in [13], will be given in brief followed by our proposed method which is to be applied for kinked boundaries.  In Section 4 we shall present a sample problem to demonstrate that the proposed method is able to produce finite element solution for general unbounded domain with tractions or prescribed values at boundaries exhibiting decay condition. We shall conclude our discussions in Section 5.
-The layout of this part of the paper is as follows; in the next section we shall explain the test procedure including the generalized formulation given in the first part of the paper.  Some sample problems are added to the section in order to give insight to operations needed especially for the new part of the formulation.  In Section 3, a brief overview of the error estimators used are given.  In Section 4 the results obtained from application of the test to heat and elasticity problem are presented and discussed.  Finally, in Section 5 we shall summarize the conclusions made throughout the paper.
+=2. PRELIMINARIES=
-In this part of the paper we shall frequently refer to relations and equation numbers in the previous part. In order to distinguish the equation numbers from the current part and the previous one, we use symbol “[[Image:draft_Samper_725107104-picture-x0000_i1025.svg|19px]] ” in the quotation of the numbers when the equation belongs to the first part, e.g. (30-[[Image:draft_Samper_725107104-picture-x0000_i1026.svg|19px]] ) refers to equation number 30 of the first part.
+In this section we prepare the preliminary tools for study of asymptotic behaviors of finite element solutions as well as error estimation schemes.
-==2. The test procedure==
+The problems of interest in this study are elliptic ones, e.g. heat transfer or elasticity problems, with generic differential equation/equations as
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math display="inline">\boldsymbol{S}^T\boldsymbol{DSu}+\boldsymbol{b}=\boldsymbol{0}</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (1-a)
+|}
+to be solved on  <math display="inline">\Omega </math> with following boundary condition
-In the first part of the paper [1], preliminary concepts related to the model problem and error definition were discussed.  The test procedure for patch locations at interior parts as well as at flat boundaries originally proposed by Babuška and co workers [9-11], was revisited. A generalization for such a patch test was also introduced in the previous part.  Here we shall give some details for the procedure of test near kinked boundaries.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math display="inline">\boldsymbol{u}=\boldsymbol{u}_0</math> on     <math display="inline">{\Gamma }_u</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (1-b)
+|}
-As is seen in the first part, the test procedure is consisting of two main parts; first part is evaluation of the asymptotic finite element solution for a given exact solution as monomials of a [[Image:draft_Samper_725107104-picture-x0000_i1027.svg|45px]] polynomial and the second part is estimation of the error and calculation of the upper and lower bounds of the effectivity index.
+and
-The asymptotic finite element solution is considered as summation of two sets of solutions, one as an asymptotic solution for interior parts of the domain and another as a boundary layer solution which plays the role of a correction to the first set when the patch used is near a flat or a kinked boundary.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math display="inline">\boldsymbol{\tilde{n}DSu}-\boldsymbol{t}=\boldsymbol{0}</math> on     <math display="inline">{\Gamma }_t</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (1-c)
+|}
-The boundary layer solution is in fact a spurious mode which decays towards the interior parts of the domain.  The asymptotic boundary layer solution near flat boundaries may be obtained using the method proposed by Babuška ''et al'' in [20] through which an eigenvalue problem is defined and solved in order to impose the decay condition in direction of normal to the boundary.  However, one may invariably use the extended formulation given in the first part by fixing one of the proportionality factors to unity.  For the case of kinked boundaries the solution is performed in a spectral form using an eigenvalue solution for each component.   It is shown that the procedure is able to produce a finite element solution for an unbounded domain [1].
-Having found the asymptotic finite element solution for a monomial, the test procedure is further continued towards evaluation of the error.  In this paper we shall examine the performance of some recovery based error estimators.  In such an error estimator new fields of gradients are constructed and the error is evaluated approximately as
+Where  <math display="inline">{\Gamma }_u\cap {\Gamma }_t=\varnothing </math> and   <math display="inline">{\Gamma }_u\cup {\Gamma }_t=\partial \Omega </math> .  In above  <math display="inline">\boldsymbol{S}</math> is a suitable operator for defining the gradients of  <math display="inline">\boldsymbol{u}</math> ,  <math display="inline">\boldsymbol{D}</math> is a suitable matrix for modulus of material and  <math display="inline">\boldsymbol{t}</math> is a traction vector associated with the problem type.  Also  <math display="inline">\boldsymbol{\tilde{n}}</math> is a matrix containing components of unit normal to the boundary arranged in an appropriate form for heat or elasticity problems.
-[[Image:draft_Samper_725107104-picture-x0000_i1028.svg|77px]] ,         [[Image:draft_Samper_725107104-picture-x0000_i1029.svg|75px]] or         [[Image:draft_Samper_725107104-picture-x0000_i1030.svg|83px]] (1)
+Discretization of the domain to a series of elements and approximation of  <math display="inline">\boldsymbol{u}</math> as  <math display="inline">\boldsymbol{u}\cong \boldsymbol{u}_h=\boldsymbol{N}{\boldsymbol{\overline{u}}}_h</math> and application of Galerkin type of weighted residual formulation in weak form (or equivalent minimization of a suitable functional) leads to following well known matrix equation
-Where parameters with star represent the enhanced-recovered solutions and those with subscript [[Image:draft_Samper_725107104-picture-x0000_i1031.svg|13px]] denote the corresponding values from asymptotic finite element solution.   Therefore the error norms required for definition of effectivity index ( 6-[[Image:draft_Samper_725107104-picture-x0000_i1032.svg|19px]] ) are written as (7-[[Image:draft_Samper_725107104-picture-x0000_i1033.svg|19px]] ).
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\int }_{\mbox{ }\Omega }\boldsymbol{B}^T\boldsymbol{DB}{\boldsymbol{\overline{u}}}_h\mbox{ }d\Omega -</math><math>{\int }_{\mbox{ }\Omega }\boldsymbol{N}^T\boldsymbol{b}\mbox{ }d\Omega -</math><math>{\int }_{\mbox{ }{\Gamma }_t}\boldsymbol{N}^T\boldsymbol{t}\mbox{ }d\Gamma =</math><math>\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (2)
+|}
-Of course the approximate error norm may be compared to the norm of the exact errors obtained from the difference of the exact gradients, pertinent to the monomial used as the exact solution, and the asymptotic finite element solution.  Considering all possible error norms from the monomials one can calculate the upper and lower bounds of the effetivity index as (19-[[Image:draft_Samper_725107104-picture-x0000_i1034.svg|19px]] ) to (22-[[Image:draft_Samper_725107104-picture-x0000_i1035.svg|19px]] ) and the robustness index as (8-[[Image:draft_Samper_725107104-picture-x0000_i1036.svg|19px]] ).
-In summary the test procedure near intersecting boundaries or even a flat boundary may be described as below
+Where  <math display="inline">\boldsymbol{N}</math> is the set of shape functions used and also  <math display="inline">\boldsymbol{B}\equiv \boldsymbol{SN}</math> .  It is noteworthy to remember that the following relation always exists between the exact solution and the finite element one
-:(1)  Select a repeatable pattern for the patch.
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math display="inline">{\int }_{\mbox{ }\Omega }\boldsymbol{B}^T\boldsymbol{DB}{\boldsymbol{\overline{u}}}_h\mbox{ }d\Omega =</math><math>{\int }_{\mbox{ }\Omega }\boldsymbol{B}^T\boldsymbol{DS}\boldsymbol{u}_{ex}\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (3)
+|}
-:(2)  Select a proper monomial with one order higher than the order of the shape functions.
-:(3)  Solve Equation (23-[[Image:draft_Samper_725107104-picture-x0000_i1037.svg|19px]] ) for [[Image:draft_Samper_725107104-picture-x0000_i1038.svg|21px]] with periodic boundary conditions (24-[[Image:draft_Samper_725107104-picture-x0000_i1039.svg|19px]] ) and minimum supports for the interior patch.
+Generally there will be a discrepancy between the exact solution, which is not available for most problems, and the finite element solution.  This deviation constitutes the error at each point
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math display="inline">\boldsymbol{e}_u=\boldsymbol{u}_{ex}-\boldsymbol{u}_h</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (4-a)
+|}
+The error distribution can also be expressed in terms of the gradients
-:(4)  Use Equation (16-[[Image:draft_Samper_725107104-picture-x0000_i1040.svg|19px]] ) to find [[Image:draft_Samper_725107104-picture-x0000_i1041.svg|16px]] for interior patches.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math display="inline">\boldsymbol{e}_{\epsilon }={\epsilon }_{ex}-{\epsilon }_h</math> or          <math display="inline">\boldsymbol{e}_{\sigma }={\sigma }_{ex}-{\sigma }_h</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (4-b)
+|}
-:(5)  Consider a [[Image:draft_Samper_725107104-picture-x0000_i1042.svg|33px]] (or [[Image:draft_Samper_725107104-picture-x0000_i1043.svg|32px]] ) macro patch. If the test procedure is to be performed near a corner ,with three intersecting flat boundaries, use a [[Image:draft_Samper_725107104-picture-x0000_i1044.svg|53px]] (or [[Image:draft_Samper_725107104-picture-x0000_i1045.svg|51px]] ) macro patch.
-:(6)  Construct the stiffness matrix of each patch in the macro patch and condense out those degrees of freedoms which are not shared with adjacent patches. Assemble the stiffness matrices to construct the macro patch stiffness matrix.
+Where  <math display="inline">{\sigma }_{ex}=\boldsymbol{D}{\epsilon }_{ex}=\boldsymbol{DS}\boldsymbol{u}_{ex}</math> and <math display="inline">{\sigma }_h=\boldsymbol{D}{\epsilon }_h=\boldsymbol{DB}{\boldsymbol{\overline{u}}}_h</math> .
-:(7)  Write the proportionality relations like (52-[[Image:draft_Samper_725107104-picture-x0000_i1046.svg|19px]] ) between the nodal values on shared faces of the patches.  For a corner patch, write the relations in all directions.  For an edge patch, however, a periodicity condition with unit proportionality factor must be used along the edge.
+If the shape functions are chosen properly, reduction of the maximum mesh size must result to reduction of errors, i.e.  <math display="inline">\boldsymbol{e}\rightarrow \boldsymbol{0}</math> as  <math display="inline">h_{max}\rightarrow 0</math> .  This is usually referred to as completeness of the set of shape functions.
-:(8)  Select a number of patches along the axes with non periodic conditions (along axes with periodic conditions just one patch is enough).
+It is a priori understood that in absence of any singularity in the domain, the ideal rate for convergence of the finite element errors is of order  <math display="inline">O(h^{P+1})</math> for  <math display="inline">\boldsymbol{e}_u</math> and  <math display="inline">O(h^P)</math> for  <math display="inline">\boldsymbol{e}_{\epsilon }</math> and  <math display="inline">\boldsymbol{e}_{\sigma }</math> where  <math display="inline">h</math> is an average mesh size for the solution.  This simple information constitutes a basis for study of error behavior and evaluation of error estimates in asymptotic form.
-:(9)  Choose a numerical integration scheme and select a series of integration points.
+.1  ''Error estimation''
-:(10)  For each integration point construct the transformation matrix as used in (53-[[Image:draft_Samper_725107104-picture-x0000_i1047.svg|19px]] ) to relate the dependent nodal values, in the macro patch,  to independent values in the master patch.
+Approximate evaluation of error can be performed through estimation of its absolute value, irrespective of the signs, or calculation of the difference between a set of enhanced answers and those of finite element method.  Most of residual based error estimators fall within the former approach and almost all recovery based error estimators employ the latter approach.
-:(11)  Evaluate [[Image:draft_Samper_725107104-picture-x0000_i1048.svg|21px]] as Equation (56-[[Image:draft_Samper_725107104-picture-x0000_i1049.svg|19px]] ).
+In an error analysis a suitable norm is usually defined to give an indication for the exact error level over the whole domain, e.g.  <math display="inline">L_2</math> norm of errors
-:(12)  For each integration point find the roots of the determinant of [[Image:draft_Samper_725107104-picture-x0000_i1050.svg|21px]] and the associative null space vectors. Choose the roots that satisfy decay condition.
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math display="inline">{\Vert \boldsymbol{e}_u\Vert }_{L_2,{\Omega }^{{'}}}=</math><math>{\left({\int }_{\mbox{ }{\Omega }^{{'}}}{\left(\boldsymbol{u}_{ex}-\boldsymbol{u}_h\right)}^T(\boldsymbol{u}_{ex}-\boldsymbol{u}_h)\mbox{ }d\Omega \right)}^{\frac{1}{2}}</math> or       <math display="inline">{\Vert \boldsymbol{e}_{\sigma }\Vert }_{L_2,{\Omega }^{{'}}}=</math><math>{\left({\int }_{\mbox{ }{\Omega }^{{'}}}{\left({\sigma }_{ex}-{\sigma }_h\right)}^T({\sigma }_{ex}-{\sigma }_h)\mbox{ }d\Omega \right)}^{\frac{1}{2}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (5)
+|}
-:(13)  Construct the transformation matrix[[Image:draft_Samper_725107104-picture-x0000_i1051.svg|19px]] in (66-[[Image:draft_Samper_725107104-picture-x0000_i1052.svg|19px]] ) or [[Image:draft_Samper_725107104-picture-x0000_i1053.svg|19px]] (76-[[Image:draft_Samper_725107104-picture-x0000_i1054.svg|19px]] ) for application of Dirichlet or Neumann conditions, respectively, to generate vectors [[Image:draft_Samper_725107104-picture-x0000_i1055.svg|13px]] or [[Image:draft_Samper_725107104-picture-x0000_i1056.svg|17px]] .
-:(14)   Construct the transformation matrix [[Image:draft_Samper_725107104-picture-x0000_i1057.svg|19px]] in (64-[[Image:draft_Samper_725107104-picture-x0000_i1058.svg|19px]] ) to generate vector of nodal values for the domain of interest, i.e. [[Image:draft_Samper_725107104-picture-x0000_i1059.svg|17px]] in (64-[[Image:draft_Samper_725107104-picture-x0000_i1060.svg|19px]] ).
+Here  <math display="inline">{\Omega }^{{'}}</math> represents a partition of the domain where the error analysis is to be performed. When evaluation of an error estimator is required, benchmark problems with exact solution are selected and the true error and the approximate one are compared.  Having evaluated approximate errors, i.e.  <math display="inline">\boldsymbol{e}_u^\ast </math> ,  <math display="inline">\boldsymbol{e}_{\epsilon }^\ast </math> or  <math display="inline">\boldsymbol{e}_{\sigma }^\ast </math> , a dimensionless parameter as the ratio of the error values, called “''effectivity index''”, is defined to show the quality of error estimates
-:(15)  Find [[Image:draft_Samper_725107104-picture-x0000_i1061.svg|28px]] and [[Image:draft_Samper_725107104-picture-x0000_i1062.svg|32px]] (or  [[Image:draft_Samper_725107104-picture-x0000_i1063.svg|48px]] if Neumann conditions are to be applied).
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math display="inline">{\theta }_{\sigma ,L_2,{\Omega }^{{'}}}=\frac{{\Vert \boldsymbol{e}_{\sigma }^\ast \Vert }_{L_2,{\Omega }^{{'}}}}{{\Vert \boldsymbol{e}_{\sigma }\Vert }_{L_2,{\Omega }^{{'}}}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (6)
+|}
-:(16)   Repeat from (j) to complete the integration procedures for two matrices in (o).
-:(17)  Evaluate [[Image:draft_Samper_725107104-picture-x0000_i1064.svg|17px]] as Equation (71-[[Image:draft_Samper_725107104-picture-x0000_i1065.svg|19px]] ).
+where
-:(18)  Apply the boundary conditions
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math display="inline">{\Vert \boldsymbol{e}_{\sigma }^\ast \Vert }_{L_2,{\Omega }^{{'}}}=</math><math>{\left({\int }_{\mbox{ }\Omega }{\left(\boldsymbol{e}_{\sigma }^\ast \right)}^T(\boldsymbol{e}_{\sigma }^\ast )\mbox{ }d\Omega \right)}^{\frac{1}{2}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (7)
+|}
-''For Dirichlet boundary conditions; '' Directly use Equation (73-[[Image:draft_Samper_725107104-picture-x0000_i1066.svg|19px]] ) to find [[Image:draft_Samper_725107104-picture-x0000_i1067.svg|24px]] for the whole domain.
-'' ''
+Obviously when  <math display="inline">\theta </math> is close to unity, the performance of the estimator is considered satisfactory.  It is expected that the effectivity index tend to unity when the mesh size is reduced sufficiently, i.e.  <math display="inline">\theta \rightarrow 1</math> as  <math display="inline">h\rightarrow 0</math> .  If such a property exists, the error estimator is called “asymptotically exact” [14]. It has been shown that only few estimators exhibit such a property (the reader may consult to references 3,4 and 10).
-''For Neumann boundary conditions; '' Evaluate [[Image:draft_Samper_725107104-picture-x0000_i1068.svg|33px]] defined in (74-[[Image:draft_Samper_725107104-picture-x0000_i1069.svg|19px]] ) and (75-[[Image:draft_Samper_725107104-picture-x0000_i1070.svg|19px]] ) and use Equation (84-[[Image:draft_Samper_725107104-picture-x0000_i1071.svg|19px]] ) or (88-[[Image:draft_Samper_725107104-picture-x0000_i1072.svg|19px]] ) to find [[Image:draft_Samper_725107104-picture-x0000_i1073.svg|24px]] for the whole domain. If Equation (84-[[Image:draft_Samper_725107104-picture-x0000_i1074.svg|19px]] ) is to be used, [[Image:draft_Samper_725107104-picture-x0000_i1075.svg|21px]] must be evaluated through Equation (82-[[Image:draft_Samper_725107104-picture-x0000_i1076.svg|19px]] ) in advance.  The prescribed nodal forces may also be calculated as the right hand side of  Equation (36-[[Image:draft_Samper_725107104-picture-x0000_i1077.svg|19px]] ).
+However, in study of asymptotic behavior of error estimators, it happens that the asymptotic effectivity indices take values close to unity.  When comparison of different types of estimators is aimed, an error estimator is said to be more robust if the effectivity index calculated is closer to unity.
-:(19)  If necessary repeat from step number (h) by increasing the number patches along the axes to ensure that no significant change happens and reasonable convergence is achieved.
+=3. ROBUSTNESS PATCH TEST=
-:(20)  Find the FEM solution though Equation (47-[[Image:draft_Samper_725107104-picture-x0000_i1078.svg|19px]] ).
+In this section we describe the test introduced by Babuška ''et al'' in 1994.  Although the procedure can be found, with many proofs for theorems and lemmas, in references [8] and [9] we repeat the core of the formulation here since it is needed for generalization of the formulation to wider set of cases.
-:(21)  Perform the recovery procedure using FEM solution from previous step.
+The test is basically a patch test since all computations at last is performed on a small group of elements with specific boundary conditions.  The whole procedure resembles to the patch tests usually employed for examining completeness of shape functions in element technology.  The reader can find the basis and history of such patch tests in reference [15].   But before explaining the robustness test, it is worthy to recall the basic features of the patch test in element technology in order to comprehend similarities and differences with the robustness patch test for error estimators.
-:(22)  Repeat from (b) to (u) for all possible monomials.
+It is well understood that if in a simple patch test an element can reproduce an exact field with order similar to the order of the shape functions then the asymptotic rate of convergence of the solution is guaranteed.  The test is in fact a finite element solution on a group of elements with boundary conditions extracted from the exact field.  The choice of boundary condition for such a test is trivial to determine and this is usually accomplished by considering prescribed values for main function, i.e. displacements, or its derivatives as prescribed tractions.
-:(23)  Choose a patch in a desired location and evaluate [[Image:draft_Samper_725107104-picture-x0000_i1079.svg|20px]] and [[Image:draft_Samper_725107104-picture-x0000_i1080.svg|16px]] (see Eqn. (17-[[Image:draft_Samper_725107104-picture-x0000_i1081.svg|19px]] ) and (18-[[Image:draft_Samper_725107104-picture-x0000_i1082.svg|19px]] )).
+It may be thought that similar strategy is applicable when asymptotic convergence of an error estimate, e.g. a superconvergent one, is to be studied.  It is true in some respects but the test will not be a complete one.   It is possible to construct a test on a group of elements with prescribed boundary conditions extracted from a polynomial, as an exact solution, with one order higher than that of the shape functions. In such a study the aim will be evaluation of capability of reproducing the exact errors.  But we note that the boundary conditions considered are not necessarily similar to those which exist between the patches of elements in a real finite element solution.  Therefore the results obtained from such a patch test may be considered just as a rough indication for capability of the error estimator.
-:(24)  Find eigenpairs of [[Image:draft_Samper_725107104-picture-x0000_i1083.svg|20px]] .
+The deficiency is removed in the patch test proposed by Babuška ''et al'' in an elegant manner.  The test is constructed by some basic assumptions for mesh configuration and smoothness of the solution.  The robustness of an error estimation, in a norm as (5) for instance, is defined by evaluation of upper and lower bounds,  <math display="inline">{\theta }_L</math> and  <math display="inline">{\theta }_U</math>, for effectivity index of (6).  Therefore the following index is defined
-:(25)  Evaluate [[Image:draft_Samper_725107104-picture-x0000_i1084.svg|16px]] from (20-[[Image:draft_Samper_725107104-picture-x0000_i1085.svg|19px]] ).
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>R=max\left(\vert 1-{\theta }_L\vert +\vert 1-{\theta }_U\vert ,\vert 1- \frac{1}{{\theta }_L}\vert +\vert 1-\frac{1}{{\theta }_U}\vert \right)</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (8)
+|}
-:(26)  Find eigenvalues of [[Image:draft_Samper_725107104-picture-x0000_i1086.svg|16px]] to determine the upper and lower bounds of the effectivity index as (22-[[Image:draft_Samper_725107104-picture-x0000_i1087.svg|19px]] ) and calculate robustness index (8-[[Image:draft_Samper_725107104-picture-x0000_i1088.svg|19px]] ).
-In following sample problem we give some numeric details.
+If the error estimator precisely reproduces the error then  <math display="inline">{\theta }_L={\theta }_U=1</math> and thus  <math display="inline">R=0</math>. However for practical use, robustness values between  <math display="inline">0.8<R<1.2</math> are considered ideal. Since the assessment is performed in an asymptotic manner, the exact solution for evaluation of the exact error is considered as a polynomial of one order higher than that of shape functions representing largest terms in a Taylor expansion of error when the size of elements, at limit, is infinitesimal.
-''Sample problem 1''
+.1 ''Assumptions''
-Supposing that a repeatable pattern of patches of elements as Figure 1 is prepared, we shall look for an asymptotic finite element solution, [[Image:draft_Samper_725107104-picture-x0000_i1089.svg|17px]] , for a scalar type of problem when a monomial as [[Image:draft_Samper_725107104-picture-x0000_i1090.svg|24px]] plays the role of exact solution and Dirichlet boundary condition is of concern. The mash pattern has been used once in the first part of this paper for solution of a general unbounded domain. Here we specifically focus on evaluation of  [[Image:draft_Samper_725107104-picture-x0000_i1091.svg|23px]] near corner of the domain since numeric details of the rest of the procedure are similar to what has been given in reference [20] for flat boundaries.
+The fundamental assumptions are as follows
-As a first step, we find [[Image:draft_Samper_725107104-picture-x0000_i1092.svg|21px]] and [[Image:draft_Samper_725107104-picture-x0000_i1093.svg|16px]] to be used in Equation (32-[[Image:draft_Samper_725107104-picture-x0000_i1094.svg|19px]] ) for boundary conditions at [[Image:draft_Samper_725107104-picture-x0000_i1095.svg|37px]] and[[Image:draft_Samper_725107104-picture-x0000_i1096.svg|37px]] .   If the nodes are numbered as Figure 1-a the nodal values of the periodic finite element solution will be
+:* It is assumed that in all directions the mesh is locally uniform and is composed of a series of repeatable pasterns, as the basic cell with dimension  <math display="inline">h</math> especially inside  <math display="inline">{\Omega }_1</math> in Figure 1. A single cell of the repeatable pattern is considered as a patch for the test.  The reader may notice that such patches have different definitions from the ones used for recovery methods.  Note that no assumption is made about the mesh outside of  <math display="inline">{\Omega }_1</math>.
-[[Image:draft_Samper_725107104-picture-x0000_i1097.svg|268px]] ,            [[Image:draft_Samper_725107104-picture-x0000_i1098.svg|279px]]
-and [[Image:draft_Samper_725107104-picture-x0000_i1099.svg|93px]] .  The independent nodal values are arranged as [[Image:draft_Samper_725107104-picture-x0000_i1100.svg|124px]] .
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="text-align: center;vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image429.png|270px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image430.png|288px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-The prescribed values of [[Image:draft_Samper_725107104-picture-x0000_i1101.svg|23px]] are calculated through equation (32-[[Image:draft_Samper_725107104-picture-x0000_i1102.svg|19px]] ). For simplicity in the step by step numerical presentation, we shall apply the boundary condition on a small domain of [[Image:draft_Samper_725107104-picture-x0000_i1103.svg|33px]] patches.
+Figure 1.  The domain, subdomain and basic cell considered for robustness patch test; (a) Three  dimensional presentation of the domain; (b) Projection of the domain on a two dimensional plane.
-Application of steps (e) to (i) leads to evaluation of [[Image:draft_Samper_725107104-picture-x0000_i1104.svg|21px]] as the one given in the sample problem of first part of the paper. Instead of evaluation of a parametric form of [[Image:draft_Samper_725107104-picture-x0000_i1105.svg|21px]] in terms of both [[Image:draft_Samper_725107104-picture-x0000_i1106.svg|19px]] and [[Image:draft_Samper_725107104-picture-x0000_i1107.svg|20px]] , one can choose a value for [[Image:draft_Samper_725107104-picture-x0000_i1108.svg|19px]] as an integration point and write [[Image:draft_Samper_725107104-picture-x0000_i1109.svg|21px]] in terms of the other parameter.  Note that the patch size in this sample problem is unity which is different from patch sizes used in [1] but this does not affect the computations pertinent to the characteristic equation and the null space.  We choose a 30 point Gauss integration rule for [[Image:draft_Samper_725107104-picture-x0000_i1110.svg|19px]] (and [[Image:draft_Samper_725107104-picture-x0000_i1111.svg|20px]] ).  For example an integration point as [[Image:draft_Samper_725107104-picture-x0000_i1112.svg|55px]] one may calculate
-[[Image:draft_Samper_725107104-picture-x0000_i1113.svg|600px]]
-from the characteristic equation (see [[Image:draft_Samper_725107104-picture-x0000_i1114.svg|64px]] in the sample problem of part I).  The null space of [[Image:draft_Samper_725107104-picture-x0000_i1115.svg|21px]] are then obtained as
+:* Study on asymptotic convergence of the solution requires that the elements sizes be close to zero and so the dimensions of patches consisting a number of such elements like patches shown in Figure 1, i.e. basic cell with dimension  <math display="inline">h</math> and the larger patches with dimensions  <math display="inline">h_0</math> and  <math display="inline">h_1</math> (<math display="inline">h_0</math> and  <math display="inline">h_1</math> are multiples of  <math display="inline">h</math>). In this regard it is also assumed that that  <math display="inline">h</math> tends to zero with a faster rate than  <math display="inline">h_0</math> (and  <math display="inline">h_1</math>). Moreover, it is assumed that the basic cells are cuboid (in three dimensions) with edges parallel to axes, though this is not necessary in general.
-'''[[Image:draft_Samper_725107104-picture-x0000_i1116.svg|120px]] ,   [[Image:draft_Samper_725107104-picture-x0000_i1117.svg|123px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1118.svg|2px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1119.svg|2px]] '''
+:* It is assumed that the exact solution is smooth at vicinity of the patch, i.e. inside macro patch  <math display="inline">{\Omega }_1</math>.  This means that the exact solution has bounded  <math display="inline">\left(p+2\right)</math> derivatives at the location of the patch (or in other words the patch is far from singular points).  Besides, not all of the derivatives with order of  <math display="inline">\left(p+1\right)</math> vanish at the location of the patch.
-We choose vectors associated with [[Image:draft_Samper_725107104-picture-x0000_i1120.svg|47px]] , i.e.[[Image:draft_Samper_725107104-picture-x0000_i1121.svg|17px]] , [[Image:draft_Samper_725107104-picture-x0000_i1122.svg|19px]] and [[Image:draft_Samper_725107104-picture-x0000_i1123.svg|20px]] , reflecting decay condition.   For each case we define the transformation matrix [[Image:draft_Samper_725107104-picture-x0000_i1124.svg|19px]] in order to construct [[Image:draft_Samper_725107104-picture-x0000_i1125.svg|13px]] , for instance for two patches along each axis (see Figure 1-b)
+:* The mesh is sufficiently refined outside  <math display="inline">{\Omega }_1</math> such that the error converges almost with the optimal rate in  <math display="inline">L_2</math> norm in  <math display="inline">{\Omega }_1</math> , i.e.:
-[[Image:draft_Samper_725107104-picture-x0000_i1126.svg|600px]]
+<math display="inline">{\Vert \boldsymbol{e}_u\Vert }_{L_2,{\Omega }_1}\leq Ch^{p+1-\lambda }</math>
-Therefore
+where  <math display="inline">C</math> is a constant independent of  <math display="inline">h</math>, and  <math display="inline">\lambda </math> representing the power of singularity, is sufficiently small. This means that the pollution error is negligible. This assumption practically holds for meshes created by adaptive algorithms which try to distribute, equally, the energy norm of the error over the elements.
-{|
+.2 ''Basic features of the test''
+Before focusing on the error, asymptotic finite element solutions should be obtained.  This means that the problem of (1) should be solved either irrespective of the main boundary conditions, when the patch is located inside, or just with condition specified when the patch located close to a boundary but still irrespective of other boundary conditions.  The information available is that the exact solution, in asymptotic sense, is of a polynomial form. If such a finite element solution is obtained the results are generalized to many problems with similar type and different conditions at boundaries.  Error evaluation comes after determining the finite element solution.
+From the smoothness assumptions it is concluded that the exact solution may be written in the form of a Taylor series around a point inside the macro-patch  <math display="inline">{\Omega }_0</math> as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1127.svg|120px]]
+| <math>u_{ex}=a_0+a_1\Delta x+a_2\Delta y+a_3\Delta z+\cdots +</math><math>a_{10}\Delta x^3+a_{11}\Delta x\Delta y^2+a_{12}\Delta x\Delta y\Delta z+</math><math>\cdots </math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1128.svg|center|196px]]
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (9)
 |}
-[[Image:draft_Samper_725107104-picture-x0000_i1129.svg|195px]]
-Similar steps are repeated for evaluation of [[Image:draft_Samper_725107104-picture-x0000_i1130.svg|19px]] in terms of [[Image:draft_Samper_725107104-picture-x0000_i1131.svg|20px]] . We evaluate matrix [[Image:draft_Samper_725107104-picture-x0000_i1132.svg|17px]] as (71-[[Image:draft_Samper_725107104-picture-x0000_i1133.svg|19px]] ) the details of which can not be given here because of numerous operations involved
-[[Image:draft_Samper_725107104-picture-x0000_i1134.svg|600px]]
+with  <math display="inline">a_i</math> being appropriate coefficients and  <math display="inline">\Delta x=x-x_0</math> and so forth for  <math display="inline">\Delta y</math> and  <math display="inline">\Delta z</math> where  <math display="inline">x_0</math>, <math display="inline">y_0</math> and  <math display="inline">z_0</math> are coordinates of a point inside the patch.   With the knowledge of asymptotic convergence of the finite element solution it follows that it is just necessary to consider a complete polynomial of order  <math display="inline">\left(p+1\right)</math> as the asymptotic exact solution.  For scalar problems, for instance
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>u_{ex}=\boldsymbol{p}^T\boldsymbol{a}</math> ,           <math>\boldsymbol{p}=\langle X^{p+1},X^PY,\cdots X^lY^mZ^n,\cdots \rangle </math> , <math display="inline">\boldsymbol{a}={\langle a_1,a_2,\cdots \cdots \rangle }^T</math> ,      <math>l+m+n=p+1</math> .
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (10)
+|}
-The solution for [[Image:draft_Samper_725107104-picture-x0000_i1135.svg|23px]] is obtained through Equation (73-[[Image:draft_Samper_725107104-picture-x0000_i1136.svg|19px]] ) using following prescribed values
-[[Image:draft_Samper_725107104-picture-x0000_i1137.svg|600px]]
+Here  <math display="inline">X</math>, <math display="inline">Y</math> and  <math display="inline">Z</math> are suitable normalized coordinates.  For testing the estimates on linear elements in three dimensional elasticity problems, for example, we write the polynomial vector as
-For example, if the nodal values for the corner patch is of concern, one may first define [[Image:draft_Samper_725107104-picture-x0000_i1138.svg|19px]] as in (64-[[Image:draft_Samper_725107104-picture-x0000_i1139.svg|19px]] )  like the following one for {|
+{| class="formulaSCP" style="width: 100%; text-align: center;"
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1140.svg|164px]]
+|
-| [[Image:draft_Samper_725107104-picture-x0000_i1141.svg|center|600px]]
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{u}_{ex}=\sum_{i=1}^3(\boldsymbol{p}^T\boldsymbol{a}_i)\boldsymbol{1}_i</math> ,              <math>\boldsymbol{p}=\langle X^2,Y^2,Z^2,XY,YZ,XZ\rangle </math> ,  <math>\boldsymbol{a}_i={\langle a_{6i-5},a_{6i-4},a_{6i-3},a_{6i-2},a_{6i-1},a_{6i}\rangle }^T</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (11)
 |}
-Thus for three pairs at point [[Image:draft_Samper_725107104-picture-x0000_i1142.svg|55px]] , vectors [[Image:draft_Samper_725107104-picture-x0000_i1143.svg|20px]] are evaluated as
+where   <math display="inline">\boldsymbol{1}_1={\langle \begin{array}{ccc}
+& 0 & 0
+\end{array}\rangle }^T</math> ,  <math display="inline">\boldsymbol{1}_2={\langle \begin{array}{ccc}
+& 1 & 0
+\end{array}\rangle }^T</math> and  <math display="inline">\boldsymbol{1}_3={\langle \begin{array}{ccc}
+& 0 & 1
+\end{array}\rangle }^T</math> .
+Now with the aid of the following trivial relation
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{u}_{ex}=\boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}+</math><math>\left(\boldsymbol{u}_{ex}-\boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (12)
+|}
+in which  <math display="inline">{\boldsymbol{\overline{u}}}_{ex}</math> is a vector of exact values at nodal points of the mesh, one may use Equation (3) and write
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }\mbox{ }\Omega }\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }d\Omega \right)\left({\boldsymbol{\overline{u}}}_h- {\boldsymbol{\overline{u}}}_{ex}\right)= </math><math> {\int }_{\mbox{ }\mbox{ }\Omega }\boldsymbol{B}^T\boldsymbol{DS}\left(\boldsymbol{u}_{ex}- \boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (13)
+|}
+Wherever the mesh is repeated, the right hand side of (13) is a periodic function since  <math display="inline">\boldsymbol{u}_{ex}-\boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}</math> is periodic.  The proof is very easy and can be accomplished by writing the monomials of (10) or (11) locally in each element and substituting in expression of  <math display="inline">\boldsymbol{u}_{ex}-\boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}</math>.  It follows that the new function is not dependent on the position of the origin of the local coordinates and therefore if the element pattern is repeated the function will be repeated as well (see [16] for more explanations).
+The repeatability property of the right hand side of Equation (13) plays a key role in reducing the size of the problem to be solved.   The periodic nature of the right hand side of (13) leads to repeatability of the left hand side of the equation.  Therefore one may write the equation on a smallest repeatable unit of the mesh, or on a larger repeatable portion consisting of the units, here known as  <math display="inline">{\Omega }^{{'}}</math>
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_h^p=</math><math>{\int }_{\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DS}\left(\boldsymbol{u}_{ex}-\right. </math><math>\left. \boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)\mbox{ }\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (14)
+|}
+in which  <math>{\boldsymbol{\overline{u}}}_h^p={\boldsymbol{\overline{u}}}_h-</math><math>{\boldsymbol{\overline{u}}}_{ex}</math>.  Equation (14) should of course be solved with a set of periodic boundary conditions.  In the next subsection we shall see that depending on the patch location  <math display="inline">{\Omega }^{{'}}</math> may be interpreted as a patch volume itself or volume of a column of patches normal to the boundary.
+It can be seen that Equation (14) is insensitive to constant values for  <math display="inline">{\boldsymbol{\overline{u}}}_h^p</math>, i.e. rigid body motion of the solution.  This implies that after finding  <math display="inline">{\boldsymbol{\overline{u}}}_h^p</math>, as will be described later, the asymptotic finite element solution is obtained as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_h={\boldsymbol{\overline{u}}}_{ex}+</math><math>{\boldsymbol{\overline{u}}}_h^p+\boldsymbol{C}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (15)
+|}
+Where  <math display="inline">\boldsymbol{C}</math> is a vector containing constant values and is evaluated from equivalence of mean values for the exact and finite element solutions [9, 10, 13, 16]
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\int }_{\mbox{ }{\Omega }^{{'}}}(\boldsymbol{N}{\boldsymbol{\overline{u}}}_h^p+</math><math>\boldsymbol{C})\mbox{ }\mbox{ }d\Omega ={\int }_{\mbox{ }{\Omega }^{{'}}}(\boldsymbol{u}_{ex}-</math><math>\boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex})\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (16)
+|}
+Having found the asymptotic finite element solution from  <math display="inline">\boldsymbol{u}_h=\boldsymbol{N}{\boldsymbol{\overline{u}}}_h</math>, the exact norm of errors and that of estimated ones are calculated.  Since the exact solution is expressed in terms of monomials with unknown coefficients  <math display="inline">\boldsymbol{a}</math>, as Equation (10) or (11), the norms of errors in squared form can be written as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\Vert \boldsymbol{e}_{\boldsymbol{\sigma }}\Vert }_{L_2,{\Omega }^{{'}}}^2=</math><math>{\int }_{\mbox{ }{\Omega }^{{'}}}{\left({\boldsymbol{\sigma }}_{ex}-{\boldsymbol{\sigma }}_h\right)}^T\left({\boldsymbol{\sigma }}_{ex}-\right. </math><math>\left. {\boldsymbol{\sigma }}_h\right)\mbox{ }d\Omega =</math><math>\boldsymbol{a}^T\boldsymbol{E}_0\boldsymbol{a}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (17)
+|}
+With analogy, the approximate error may also be written as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\Vert \boldsymbol{e}_{\boldsymbol{\sigma }}^\ast \Vert }_{L_2,{\Omega }^{{'}}}^2=</math><math>\boldsymbol{a}^T\boldsymbol{Ea}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (18)
+|}
+This makes it easy to find the upper and lower bounds of the squared effectivity index
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\theta }_{\boldsymbol{\sigma },L_2,{\Omega }^{{'}}}^2=</math><math>\frac{\boldsymbol{a}^T\boldsymbol{Ea}}{\boldsymbol{a}^T\boldsymbol{E}_0\boldsymbol{a}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (19)
+|}
+by evaluation of eigenvalues of the following matrix
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{\tilde{E}}=\boldsymbol{H}^{-T}\boldsymbol{E}\boldsymbol{H}^{-1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (20)
+|}
+Where the matrix  <math display="inline">\boldsymbol{H}</math> is defined as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{H}={\left(\boldsymbol{E}_0\right)}^{\frac{1}{2}}=</math><math>\boldsymbol{\Psi }{\boldsymbol{\Lambda }}^{\frac{1}{2}}{\boldsymbol{\Psi }}^T</math> ,               <math>\boldsymbol{E}_0=\boldsymbol{H}^T\boldsymbol{H}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (21)
+|}
+In above equation  <math display="inline">\Psi </math> and  <math display="inline">\Lambda </math> are normalized eigenvectors and diagonal matrix of eigenvalues of  <math display="inline">\boldsymbol{E}_0</math>, respectively.  Hence
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\left({\lambda }_{min}\right)}^{\frac{1}{2}}<\theta <{\left({\lambda }_{max}\right)}^{\frac{1}{2}}</math> , <math>{\theta }_L={\left({\lambda }_{min}\right)}^{\frac{1}{2}}</math> ,  <math>{\theta }_U={\left({\lambda }_{max}\right)}^{\frac{1}{2}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (22)
+|}
+It should be noted that relations (12) to (16) are useful for reducing the size of the problem only when a sort of periodicity exists in the solution.  This will be the case as long as the exact solution is expressed in terms of monomials of order  <math display="inline">p+1</math> and the mesh pattern is repeatable.  When the patch is far enough from the main boundary conditions, i.e. located at interior of the domain, the periodicity may be assumed for all directions. This allows us to reduce the size of the problem to a domain as small as a patch.  In the case that patch is located near a boundary, the repeatability of the mesh is lost in normal direction to the boundary and thus the solution loses its periodicity in that direction.  This latter case needs additional treatment which we shall briefly explain in the forthcoming subsections.  In the case that the patch is located at a corner of the main domain, the periodicity will be lost in all directions.  We shall also refer to this case after revisiting the patch test near a single flat boundary.
+''3.3 Patch test at interior parts of the domain''
+As mentioned above, when patch is located at interior parts of the domain, the patch pattern is considered repeatable in all directions.  Now since the source term in the right hand side of (14) is repeatable, the finite element solution will be periodic.   This allows us to write Equation (14) just on a single patch as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }{\Omega }_p}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_h^p=</math><math>{\int }_{\mbox{ }{\Omega }_p}\boldsymbol{B}^T\boldsymbol{DS}\left(\boldsymbol{u}_{ex}- \boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (23)
+|}
+which is to be solved with periodic boundary conditions i.e.
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+|
 {|
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1144.svg|132px]]
+| <math>{\boldsymbol{\overline{u}}}_{h(X+L_x)}^p={\boldsymbol{\overline{u}}}_{h(X)}^p</math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1145.svg|center|197px]]
+| <math>{\boldsymbol{\overline{u}}}_{h(Y+L_y)}^p={\boldsymbol{\overline{u}}}_{h(Y)}^p</math>
+|}
+<math>{\boldsymbol{\overline{u}}}_{h(Z+L_z)}^p={\boldsymbol{\overline{u}}}_{h(Z)}^p</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (24)
 |}
-[[Image:draft_Samper_725107104-picture-x0000_i1146.svg|196px]]
-Evaluation of [[Image:draft_Samper_725107104-picture-x0000_i1147.svg|20px]] for other integration points and use of [[Image:draft_Samper_725107104-picture-x0000_i1148.svg|17px]] and [[Image:draft_Samper_725107104-picture-x0000_i1149.svg|28px]] as given above in (73-[[Image:draft_Samper_725107104-picture-x0000_i1150.svg|19px]] ) leads to the following values for the boundary layer nodal values
-[[Image:draft_Samper_725107104-picture-x0000_i1151.svg|215px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1152.svg|163px]] ,
+with  <math display="inline">L_x</math>,  <math display="inline">L_y</math> and  <math display="inline">L_z</math> being the patch dimensions (period length) along the three directions.  Having found the periodic solutions, the asymptotic finite element solution is evaluated from (15) and (16).  The test procedure will be completed by estimating the errors and calculation of upper and lower bounds of the effectivity index as indicated in Equations (17) to (22).
-[[Image:draft_Samper_725107104-picture-x0000_i1153.svg|171px]] ,             [[Image:draft_Samper_725107104-picture-x0000_i1154.svg|117px]]
+''3.4 Patch test at parts near flat boundaries ''
-The ninth nodal value in the corner patch may be obtained from an inverse static condensation as [[Image:draft_Samper_725107104-picture-x0000_i1155.svg|117px]] .   Now summation of the nodal values as [[Image:draft_Samper_725107104-picture-x0000_i1156.svg|144px]] gives
+To extend the formulation to more general cases, here we revisit the procedure near flat boundaries. Now supposing that a flat boundary exists at  <math display="inline">z=0</math> and the patch pattern is repeatable in all directions. Equation (14) may be rewritten as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_{hB}^p=</math><math>{\int }_{\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DS}\left(\boldsymbol{u}_{ex}- \boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (25)
+|}
+Clearly the periodicity of solution is lost along z axis since the repeatability of the patch pattern is terminated at  <math display="inline">z=0</math> although the source term is still periodic.  Hence an additional solution is considered, as boundary layer type, to those obtained for interior parts
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{hB}^p=({\boldsymbol{\overline{u}}}_h^p+</math><math>\boldsymbol{C})+{\boldsymbol{\overline{u}}}_{h0}^p</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (26)
+|}
+Substituting (26) in (25) and noting that (14) holds for solutions pertaining to interior patches, then provided that no extra forces is imposed from boundaries, i.e. when dealing with Dirichlet boundary conditions,  the following integral equation is obtained
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_{h0}^p=\boldsymbol{0}\mbox{ }\mbox{ }\mbox{on }\mbox{ }{\Omega }^{{'}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (27)
+|}
+It is clear that as we go far from the boundary, i.e.  <math display="inline">z\rightarrow \infty </math>, we should obtain  <math display="inline">{\boldsymbol{\overline{u}}}_{hB}^p\rightarrow {\boldsymbol{\overline{u}}}_h^p+</math><math>\boldsymbol{C}</math>.  Therefore Equation (27) must be solved with two periodic conditions along  <math display="inline">x</math> and  <math display="inline">y</math> and following decay condition
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p\rightarrow 0</math> (or  <math>\boldsymbol{B}{\boldsymbol{\overline{u}}}_{h0}^p\rightarrow 0</math> )       when      <math>\frac{d_z}{L_z}\rightarrow \infty </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (28)
+|}
+Where  <math display="inline">d_z</math> denotes the distance along  <math display="inline">z</math>.  It may be noted that decay condition on  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> is stronger than the one defined for its gradients, i.e.  <math display="inline">\boldsymbol{B}{\boldsymbol{\overline{u}}}_{h0}^p</math>.  In fact application of decay condition in latter form may lead to constant values of  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> for large values of  <math display="inline">z</math> which is still consistent with formulation for interior parts since constant values are considered a priori.   This is true basically when we are dealing with asymptotic finite element solutions in which the solutions for two fictitious domains are added together as Equation (26).  However, when a general unbounded domain is to be solved, constant values for the solution at far ends of the domain may have no physical meaning and thus in such situations the first form of decay condition must be satisfied rigorously.  In a sample problem we shall give solution to an unbounded domain in order to verify the method to be described in sequel.
+The boundary conditions considered in this study are
+:1. ''Dirichlet boundary condition: '' In this type we apply
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_h={\boldsymbol{\overline{u}}}_{ex}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (29)
+|}
+Substitution of (26) instead of the last two terms in right hand of (15) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_h={\boldsymbol{\overline{u}}}_{ex}+</math><math>{\boldsymbol{\overline{u}}}_{hB}^p</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (30)
+|}
+and application of (29) leads to following condition
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{hB}^p\vert _{{\Gamma }_D}=</math><math>\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (31)
+|}
+From which it follows that
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p\vert _{{\Gamma }_D}=</math><math>-({\boldsymbol{\overline{u}}}_h^p+\boldsymbol{C})\vert _{{\Gamma }_D}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (32)
+|}
+If homogeneous Dirichlet condition is of concern then the monomials in (10) or (11) should be
+chosen accordingly.
+:2. ''Neumann boundary condition:'' In this type of condition the tractions extracted from the exact solution are applied to the boundary of the domain, i.e.
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{t}_0=\boldsymbol{\tilde{n}DS}\boldsymbol{u}_{ex}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (33)
+|}
+to be substituted in (1-c).  If Equation (2) is rewritten with  <math display="inline">{\boldsymbol{\overline{u}}}_h</math> replaced with its equivalent from (26) and (27) then
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}[{\boldsymbol{\overline{u}}}_{ex}+</math><math>({\boldsymbol{\overline{u}}}_h^p+\boldsymbol{C})+{\boldsymbol{\overline{u}}}_{h0}^p]\mbox{ }\mbox{ }d\Omega -</math><math>{\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{N}^T\boldsymbol{b}\mbox{ }d\Omega -</math><math>{\int }_{\mbox{ }\mbox{ }{\Gamma }_t}\boldsymbol{N}^T\boldsymbol{t}_0\mbox{ }d\Gamma =</math><math>\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (34)
+|}
+Where
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{b}=-\boldsymbol{S}^T\boldsymbol{DS}\boldsymbol{u}_{ex}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (35)
+|}
+It follows that
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_{h0}^p=</math><math>{\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{N}^T\boldsymbol{S}^T\boldsymbol{DS}\boldsymbol{u}_{ex}d\Omega +</math><math>{\int }_{\mbox{ }\mbox{ }{\Gamma }_t}\boldsymbol{N}^T\boldsymbol{\tilde{n}DS}\boldsymbol{u}_{ex}d\Gamma -</math><math>{\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }[{\boldsymbol{\overline{u}}}_{ex}+</math><math>{\boldsymbol{\overline{u}}}_h^p]\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (36)
+|}
+In fact we only need to consider elements near boundary since Equation (27) holds for interior parts. Therefore
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }{{\Omega }^{{'}}}_{bnd}}\boldsymbol{B}_{bnd}^T\boldsymbol{DB}\mbox{ }d\Omega \right){\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd}=\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (37)
+|}
+In above  <math display="inline">{{\Omega }^{{'}}}_{bnd}</math> denotes volume of all elements in the column having at least one node on the boundary,  <math display="inline">{\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd}</math> denotes associative subset of  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math>, and  <math display="inline">\boldsymbol{F}_{prs}</math> represents the right hand side of Equation (36).
+Now we are ready to apply the boundary condition to a problem asymptotically modeled with a periodic mesh pattern and a single flat boundary.   Considering a 3D mesh of patches like Figure 2-a, we assume that the boundary exists at plane  <math display="inline">z=0</math>.  The mesh pattern is considered repeatable along <math display="inline">x</math>, <math display="inline">y</math> and  <math display="inline">z</math>, with  <math display="inline">-\infty <x<\infty </math>,  <math display="inline">-\infty <y<\infty </math> and  <math display="inline">0<z<\infty </math>.  A column of patches along  <math display="inline">z</math> is selected.  For simplicity we use Figure 2-b instead of Figure 2-a as a two dimensional presentation of the three dimensional column.  Note that we are aiming at solving Equation (27) or (37) in which  <math display="inline">{\Omega }^{{'}}</math> is volume of the column.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image431.png|264px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image432.png|312px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
+Figure 2.  A column with infinite length; (a) Three dimensional presentation of the column with one end at boundary consisting of repeatable patterns as basic cells, (b)  Two dimensional presentation of the infinite column with repeatable pattern and the numbering used for different levels.
+After numbering the patch levels consecutively towards the middle of the domain as shown in Figure 2-b, in the matrix equation of each patch the degrees of freedoms not falling on interfaces of patches are statically condensed.  The matrix equation of the asymptotic finite element solution pertaining to the remaining degrees of freedoms is assembled as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\left[\begin{array}{ccccccccc}
+\boldsymbol{k}_{0,0} & \boldsymbol{k}_{0,1} &  &  &  &  &  &  & \\
+\boldsymbol{k}_{1,0} & \boldsymbol{k}_{1,1} & \boldsymbol{k}_{1,2} &  &  &  &  &  & \\
+ & \boldsymbol{k}_{2,1} & \boldsymbol{k}_{2,2} & \boldsymbol{k}_{2,3} &  &  &  &  & \\
+ &  &  &  & \ddots  &  &  &  & \\
+ &  &  &  &  & \boldsymbol{k}_{n,n-1} & \boldsymbol{k}_{n,n} & \boldsymbol{k}_{n,n+1} & \\
+ &  &  &  &  &  &  &  & \ddots
+\end{array}\right]}_{\infty \times \infty }\cdot {\left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,0}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,2}\\
+\vdots \\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n-1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n+1}\\
+\vdots
+\end{array}\right\}}_{\infty \times 1}={\left\{\begin{array}{c}
+\boldsymbol{F}_0\\
+\boldsymbol{0}\\
+\boldsymbol{0}\\
+\vdots \\
+\boldsymbol{0}\\
+\boldsymbol{0}\\
+\boldsymbol{0}\\
+\vdots
+\end{array}\right\}}_{\infty \times 1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (38)
+|}
+Where  <math display="inline">\boldsymbol{k}_{m,n}=\boldsymbol{k}_{n,m}^T</math> is the relative stiffness matrix of  <math display="inline">m</math> th and  <math display="inline">n</math> th levels, and  <math>{\boldsymbol{\overline{u}}}_{h0}^{p,k}</math> is the set of nodal values at level  <math display="inline">k</math>.
+It can be seen that each row of the matrix represents the relation between levels  <math display="inline">n-1</math>,  <math display="inline">n</math> and  <math display="inline">n+1</math> except for the first level.   Also, since the patch pattern is the same for all levels, the out of diagonal matrices, below and above the diagonal, are similar.  Extracting the  <math display="inline">n</math> th row of the matrix Equation (38), one can show it as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left[\begin{array}{ccc}
+\boldsymbol{A} & \boldsymbol{B} & \boldsymbol{C}
+\end{array}\right]\cdot \left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n-1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n+1}
+\end{array}\right\}=\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (39)
+|}
+By appropriate arrangement we have
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left[\begin{array}{cc}
+\boldsymbol{A} & \boldsymbol{B}\\
+\boldsymbol{0} & \boldsymbol{I}
+\end{array}\right]\cdot \left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n+1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}
+\end{array}\right\}=\left[\begin{array}{cc}
+\boldsymbol{0} & -\boldsymbol{C}\\
+\boldsymbol{I} & \boldsymbol{0}
+\end{array}\right]\cdot \left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n-1}
+\end{array}\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (40)
+|}
+Supposing that the coefficient matrix in the left hand side of (40) is invertible, it follows that
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n+1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}
+\end{array}\right\}=\left[\begin{array}{cc}
+-\boldsymbol{A}^{-1}\boldsymbol{B} & -\boldsymbol{A}^{-1}\boldsymbol{C}\\
+\boldsymbol{I} & \boldsymbol{0}
+\end{array}\right]\cdot \left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n-1}
+\end{array}\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (41)
+|}
+One may write a relation between two sets of such defined vectors for two levels, say levels  <math display="inline">m</math> and  <math display="inline">n</math>, as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,n+1}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,n}
+\end{array}\right\}={\left[\begin{array}{cc}
+-\boldsymbol{A}^{-1}\boldsymbol{B} & -\boldsymbol{A}^{-1}\boldsymbol{C}\\
+\boldsymbol{I} & \boldsymbol{0}
+\end{array}\right]}^{n-m}\cdot \left\{\begin{array}{c}
+{\boldsymbol{\overline{u}}}_{h0}^{p,m}\\
+{\boldsymbol{\overline{u}}}_{h0}^{p,m-1}
+\end{array}\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (42)
+|}
+Solution of above recursive relation can be obtained by transformation of the variables to a space with bases of normalized eigenvectors of the coefficient matrix in Equation (41).  Invariably, the eigenvalue problem may be solved using matrices in Equation (40), as a generalized eigenvalue problem, to avoid problems encountering for inversion of  <math display="inline">\boldsymbol{A}</math>.  Therefore following generalized eigenequation is defined
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left(\left[\begin{array}{cc}
+\boldsymbol{A} & \boldsymbol{B}\\
+\boldsymbol{0} & \boldsymbol{I}
+\end{array}\right]-\mu \left[\begin{array}{cc}
+\boldsymbol{0} & -\boldsymbol{C}\\
+\boldsymbol{I} & \boldsymbol{0}
+\end{array}\right]\right)\cdot \boldsymbol{w}=\boldsymbol{0}</math> or         <math display="inline">\left(\boldsymbol{Q}-\mu \boldsymbol{M}\right)\boldsymbol{w}=</math><math>\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (43)
+|}
+Solution of the eigenequation of (43) results in a series of  eigenvalues  <math display="inline">{\mu }_i</math> and eigenvectors  <math display="inline">\boldsymbol{w}_i</math>, as many as two times of DOFs at each level, from which those with decay property, as  <math display="inline">\vert {\mu }_i\vert \leq 1</math>, are nominated for construction of the asymptotic finite element solution.  We note that in the new space constructed by series of  <math display="inline">\boldsymbol{w}_i</math> the recursive relation for each base vector at two different levels,  <math display="inline">m</math> and  <math display="inline">n</math>, is written as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{w}_i^n={\mu }_i^{n-m}\boldsymbol{w}_i^m</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (44)
+|}
+Also each eigenvector is of the form
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{w}_i=\left\{\begin{array}{c}
+{\mu }_i{\boldsymbol{\varphi }}_i\\
+{\boldsymbol{\varphi }}_i
+\end{array}\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (45)
+|}
+Now assume that the column with infinite length is cut and the boundary conditions are to be imposed at current level with nodal values as  <math display="inline">{\varphi }_i</math>. If the levels are renumbered accordingly then for  <math display="inline">n</math> levels farther the correction to the periodic finite element solution is evaluated as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^{p,n}=\sum_{i=1}^{n_c}c_i{\mu}_i{}^n{\boldsymbol{\varphi }}_i</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (46)
+|}
+In which  <math display="inline">n_c</math> is the number of DOFs at each level and  <math display="inline">c_i</math> denotes the unknown coefficient for the mode number  <math display="inline">i</math> to be determined from the boundary condition specified in (31) or (37).
+Having found the correction required for boundaries,  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> , the periodic finite element solution may be calculated as Equation (26).  Finally, the asymptotic finite element solution is obtained as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{u}_h=\boldsymbol{N}\left({\boldsymbol{\overline{u}}}_{ex}+ {\boldsymbol{\overline{u}}}_h^p+{\boldsymbol{\overline{u}}}_{h0}^p+\right. </math><math>\left. \boldsymbol{C}\right)</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (47)
+|}
+Now one can apply the error estimation procedure to study its asymptotic convergence behavior.
+As is seen the formulation is suitable for flat boundaries only.  In case that the patch is located near a kinked boundary the periodicity of the solution is lost along to directions normal to the edge line.  The problem becomes worse when a corner patch near intersection of three planes, is considered.  In this latter case the periodicity is lost along all directions. In the next section we shall extend the formulation to such wider range of applications.
+''3.5 Generalization of the test for application to patches near kinked boundaries''
+In this section we shall give a more general formulation for evaluation of the asymptotic finite element solutions near edges or corners of a three dimensional domain.
+We assume that two intersecting planes exist at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> making a straight edge along  <math display="inline">z</math> axis (Figure 3).  A repeatable mesh pattern is again constructed.  Since the problem is studied in asymptotic form, the mesh is considered to be consisting of infinite number of patches along the three axes.  Equation (25) is again written
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_{hB}^p=</math><math>{\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DS}\left(\boldsymbol{u}_{ex}- \boldsymbol{N}{\boldsymbol{\overline{u}}}_{ex}\right)\mbox{ }d\Omega </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (48)
+|}
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 54%;"|                   [[Image:draft_Samper_725107104-image433.png|234px]]
+|  style="vertical-align: top;width: 47%;"|               [[Image:draft_Samper_725107104-image434.png|216px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
+Figure 3.  A typical case for patch location near two intersecting flat boundaries; (a) The domain shaded in gray is an unbounded domain along  <math display="inline">x</math> and  <math display="inline">y</math> axes and is of a single-cell-width  thickness, (b)  Two dimensional presentation of patch location near kinked boundaries which is similar to a corner point.
+Here  <math display="inline">{\Omega }^{{'}}</math> represents an infinite domain with thickness of a patch along  <math display="inline">z</math> axes since the periodicity of the solution retains just along  <math display="inline">z</math>.  To avoid confronting difficulties for three dimensional presentation of the problem, we look in to the problem as a two dimensional one, i.e. its projection on  <math display="inline">x</math>
+<math display="inline">y</math>
+plane (see Figure 3-b).   Again we assume that the asymptotic periodic solution is consisting of two parts as Equation (26) substitution of which in (48) gives similar equation as (27).  Clearly the solution must be found satisfying Dirichlet or Neumann boundary conditions similar to (30) or (37) and two decay conditions along  <math display="inline">x</math> and  <math display="inline">y</math>. Thus the statement of the problem is to solve
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left({\int }_{\mbox{ }\mbox{ }{\Omega }^{{'}}}\boldsymbol{B}^T\boldsymbol{DB}\mbox{ }d\Omega \right){\boldsymbol{\overline{u}}}_{h0}^p=\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (49)
+|}
+with boundary conditions (also see (32) and (37))
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p\vert _{{\Gamma }_D}=-({\boldsymbol{\overline{u}}}_h^p+\boldsymbol{C})\vert _{{\Gamma }_D}</math> or       <math>{\int }_{\mbox{ }\mbox{ }{\Omega }^'}\boldsymbol{B}^T\boldsymbol{DB}{\boldsymbol{\overline{u}}}_{h0}^pd\Omega =\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (50)
+|}
+at  <math display="inline">x=0</math> or  <math display="inline">y=0</math>, periodic conditions along  <math display="inline">z</math> and following decay conditions
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p\rightarrow 0</math> (or  <math>\boldsymbol{B}{\boldsymbol{\overline{u}}}_{h0}^p\rightarrow 0</math>)       when      <math>\frac{d_x}{L_x}\rightarrow \infty </math> and      <math>\frac{d_y}{L_y}\rightarrow \infty </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (51)
+|}
+Where here again  <math display="inline">d_x</math> and  <math display="inline">d_y</math> denote distances along  <math display="inline">x</math> and  <math display="inline">y</math>.
+The key point of the solution is isolation of a patch from  [[Image:draft_Samper_725107104-image205.png|18px]] , as before, and writing suitable relations between the DOFs. Supposing that the solution may be found spectrally in an appropriate space, some proportionality relations are assumed between adjacent patches.  The reason for assuming such relations may be explored in Appendix.  For instance in Figure 4-a following relations are considered
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+|[[Image:draft_Samper_725107104-image206.png|54px]] ,      [[Image:draft_Samper_725107104-image207.png|54px]] ,       <math display="inline">{\boldsymbol{\overline{u}}}_g={\mu }_1{\boldsymbol{\overline{u}}}_e</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (52-a)
+|}
+for  [[Image:draft_Samper_725107104-image209.png|6px]] direction and
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+|<math>{\boldsymbol{\overline{u}}}_c={\mu }_2{\boldsymbol{\overline{u}}}_h</math>,       <math>{\boldsymbol{\overline{u}}}_j={\mu }_2{\boldsymbol{\overline{u}}}_c</math>,       <math>{\boldsymbol{\overline{u}}}_l={\mu }_2{\boldsymbol{\overline{u}}}_j</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (52-b)
+|}
+for  <math display="inline">y</math> direction.  Similar relation is valid for any two nodes with distance  <math display="inline">L_x</math> or  <math display="inline">L_y</math> for example
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+|<math>{\boldsymbol{\overline{u}}}_f={\mu }_1{\boldsymbol{\overline{u}}}_d</math> ,       <math>{\boldsymbol{\overline{u}}}_k={\mu }_2{\boldsymbol{\overline{u}}}_i</math> ,      <math>{\boldsymbol{\overline{u}}}_n={\mu }_1{\boldsymbol{\overline{u}}}_m</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (52-c)
+|}
+and so forth. In above  <math display="inline">\boldsymbol{\overline{u}}</math> denotes spectral nodal values in a new space to be defined later.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="text-align: center;vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image435.png|204px]]
+|  style="text-align: center;vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image436.png|210px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
+Figure 4. A macro patch consisting of  <math display="inline">3\times 3</math> patches; (a) The proportionality relations are written as (52), between nodal values along both directions, (b) Schematic presentation of nodes in the macro patch nominated for satisfaction of the boundary condition.
+Now the solution starts with writing all relations similar to (52) for a patch and all surrounding ones, i.e. nine patches in Figure 4, in a matrix form.  It can be seen that all nodal values may be written in terms of few of them in central patch.  Note that the nodal values inside the central patch must be condensed before writing the proportionality relations.   Thus the remaining nodal values will be those located at two edges of patches, i.e. left and bottom edges.  Hence if we focus on all similar DOFs at the two similar sides of the patches, the following matrix relation containing proportionality properties, like relations (52), can be written
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_D=\boldsymbol{T}_{\left({\mu }_1,{\mu }_2\right)}{\boldsymbol{\overline{u}}}_I</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (53)
+|}
+Where  <math display="inline">{\boldsymbol{\overline{u}}}_D</math> and  <math display="inline">{\boldsymbol{\overline{u}}}_I</math> denote dependent and independent nodal values in macro patch of  <math display="inline">3\times 3</math> in Figure 4-a.  The independent nodal values in Figure 4-a may be considered as those falling in the dash line area. The reader may notice that the relation may also be written in a  [[Image:draft_Samper_725107104-image222.png|30px]] macro patch.  Here again the matrix equation obtained from (49) may be presented as (38). To avoid difficulties pertaining to numbering of the patches along two or three axes, we use an abstract notation with rearranged form as the following one
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left[\begin{array}{cc}
+\boldsymbol{K}_{II} & \boldsymbol{K}_{ID}\\
+\boldsymbol{K}_{DI} & \boldsymbol{K}_{DD} \end{array}\right]\cdot \left\{ \begin{array}{c} \boldsymbol{\overline{u}}_I\\ \boldsymbol{\overline{u}}_D   \end{array}\right\} = \left\{ \begin{matrix} \boldsymbol{0} \\ \boldsymbol{f}\end{matrix}\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (54)
+|}
+The first row of the above matrix equation represents a generic row of a matrix equation similar to (38) and the second row represents the remaining of the equations most of which have zero right hand side, i.e. components of  <math display="inline">\boldsymbol{f}</math>, except for the ones corresponding to nodes at Neumann boundaries.   Substituting (53) in (54), the first set of equations gives
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{Q^{{'}}}}_{\left({\mu }_1,{\mu }_2\right)}{\boldsymbol{\overline{u}}}_I=</math><math>\boldsymbol{0}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (55)
+|}
+In which
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{Q^{{'}}}}_{\left({\mu }_1,{\mu }_2\right)}=</math><math>\boldsymbol{K}_{II}+\boldsymbol{K}_{ID}\boldsymbol{T}_{\left({\mu }_1,{\mu }_2\right)}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (56)
+|}
+Equation (55) means that  <math display="inline"> {\boldsymbol{\overline{u}}}_I</math>   may be constructed from vectors in null space of  <math display="inline">\boldsymbol{Q^{{'}}}</math>.  To evaluate such null space the determinant of   <math display="inline">\boldsymbol{Q}^'</math> is set to be zero
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\vert {\boldsymbol{Q^{{'}}}}_{\left({\mu }_\boldsymbol{1},{\mu }_2\right)}\vert =</math><math>f({\mu }_1,{\mu }_2)=0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (57)
+|}
+It may be observed that two unknowns,  <math display="inline">{\mu }_1</math> and  <math display="inline">{\mu }_2</math>, appear in one equation.  This means that numerous pairs may be found satisfying Equation (57).  To solve the problem one may choose a value for one unknown and find the other through the characteristic equation.  Hence
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\mu }_2=g({\mu }_1)</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (58)
+|}
+To satisfy the decay condition (51) we define another function as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\mu }_2=h({\mu }_1)</math>
+| <math>h=\left\{g\mbox{ }\mbox{ }\mbox{ }\vert \mbox{ }\mbox{ }\vert g({\mu }_1)\vert \leq 1,\mbox{ }\mbox{ }\vert {\mu }_1\vert \leq 1\right\}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (59)
+|}
+Generally, several relations like (59) exist since characteristic equation (57) possesses several roots for a given  <math display="inline">{\mu }_1</math>.  Therefore
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\left({\mu }_2\right)}_i=h_i({\mu }_1)</math> when       <math>\vert h_i({\mu }_1)\vert \leq 1</math>,         <math>\vert {\mu }_1\vert \leq 1</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (60)
+|}
+Now supposing that  <math display="inline">{\phi }_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_1)}</math> belongs to the null space of  <math display="inline">\boldsymbol{Q^{{'}}}</math>, one may write the non trivial solution of (55) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_I=\int \left(\sum c_{i\mbox{ }\mbox{ }({\mu }_1)}{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_1)}\right)d{\mu }_1</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (61)
+|}
+Where integration is taken over all feasible values of  <math display="inline">{\mu }_1</math>. In Equation (61) coefficients  <math display="inline">c_{i({\mu }_1)}</math> are to be determined from boundary conditions at  <math display="inline">x=0</math> or  <math display="inline">y=0</math>.  The number of inner summation is dependent on the number of   <math display="inline">{\phi }_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_1)}</math> for which  <math display="inline">{\mu }_1</math> and  <math display="inline">{\mu }_2</math> satisfy (59) and may vary when  <math display="inline">{\mu }_1</math> falls within different regions.
+Since <math display="inline">{\mu }_1</math> may take complex values the symbolic integral of (61) is taken over a feasible area on Gaussian plane.  Because the decay condition imposes the condition of   <math display="inline">\vert {\mu }_1\vert \leq 1</math>, the feasible area is a circle with radios of unity and center at origin.  However, we note that for problems with real differential equation, the exact solution and the finite element one should be real.  This means that from values in feasible domain those leading to real answers should be selected implying that the complex values, and the associated vectors, should be in the form of conjugate pairs.  Thus the integration in (61) should in fact be taken over real values for  <math display="inline">{\mu }_1</math> while letting  <math display="inline">{\mu }_2</math> to be complex. To ensure that all real values are covered by  <math display="inline">{\mu }_1</math> and  <math display="inline">{\mu }_2</math> Equation (61) may be rewritten as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_I={\int }_{\mbox{ }0}^{\mbox{ }1}\left(\sum c_{i\mbox{ }\mbox{ }({\mu }_1)}^1{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_1)}^{\mbox{ }1}\right)\mbox{ }d{\mu }_1+</math><math>{\int }_{\mbox{ }0}^{\mbox{ }1}\left(\sum c_{i\mbox{ }\mbox{ }({\mu }_2)}^2{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_2)}^{\mbox{ }2}\right)\mbox{ }d{\mu }_2</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (62)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_I=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left(\sum c_{i\mbox{ }\mbox{ }({\mu }_k)}^k{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)\mbox{ }d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (63)
+|}
+where of course it is needed to find  <math display="inline">{\mu }_1</math> in terms of  <math display="inline">{\mu }_2</math> as Equation (58) and (59). Note that the negative values have been excluded from the integration domain because of two reasons. Firstly, change of the sign of the proportionality factor represents an oscillatory solution, i.e. oscillation within a two patch length, which may not be so desirable except for when the boundary condition has a similar oscillatory nature. Secondly, even when one factor is taken as a positive value the other, obtained as the roots of the characteristic equation, can still be negative (or conjugate pairs with negative real parts). Therefore, since both  <math display="inline">{\mu }_1</math> and <math display="inline">{\mu }_2</math> are taken as integral argument, negative values for one,  <math display="inline">{\mu }_1</math> for instance, are partially covered by the values obtained form insertion of the other factor in the equation,  <math display="inline">{\mu }_2</math> for example.  Our experience shows that this form of integration is sufficient and convergent even when an oscillatory boundary condition exists and in fact it helps to reduce the number of undesirable oscillatory modes which are to be removed from the solution.
+The above integrals may be taken numerically by a Gauss quadrature rule, for instance, excluding  <math display="inline">\mu =0</math> from the series of integration points.
+It may be noticed that  <math display="inline">\vert{\mu }_1\vert=1</math> (or  <math display="inline">\vert {\mu }_2\vert =1</math> ) does not correspond to decay condition but in fact periodic one.  But this is of no importance because of two reasons.  Firstly, if such periodic condition is not excluded from the asymptotic solution, it will not be independent of  <math display="inline">{\boldsymbol{\overline{u}}}_h^p</math> initially considered for the interior of the domain and thus the formulation will inherently include both solutions, i.e. periodic ones and the correction for the boundaries.  Secondly, if the boundary conditions specified at  <math display="inline">x=0</math> or  <math display="inline">y=0</math> inherently reflect the decay condition, i.e. vanish for large coordinates, then such a non decay solution vanishes automatically.
+Supposing that the coefficient values are available we construct the finite element solution  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> for the whole domain through defining characteristic vector  <math display="inline">\boldsymbol{w}_{i\mbox{ }\mbox{ }({\mu }_k)}^k</math> from its smallest subset, i.e.  <math display="inline">{\phi }_{\mbox{ }\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math> , using a transformation matrix like  <math display="inline">{\boldsymbol{T^{{'}}}}_{\left({\mu }_k,h({\mu }_k)\right)}^{\mbox{ }\mbox{ }k}</math> containing proportionality factors
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{w}_{i\mbox{ }\mbox{ }({\mu }_k)}^k=</math><math>{\boldsymbol{T^{{'}}}}_{\left({\mu }_k,h({\mu }_k)\right)}^{\mbox{ }\mbox{ }k}\mbox{ }{\phi }_{\mbox{ }\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (64)
+|}
+Clearly the dimension of  [[Image:draft_Samper_725107104-image258.png|12px]] depends upon the number of patches considered along two axes. Therefore correction to finite element solution due to presence of boundary condition is evaluated as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left(\sum c_{i({\mu }_k)}^k\boldsymbol{w}_{i({\mu }_k)}^k\right)d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (65)
+|}
+Now the question regarding the values of coefficients, <math display="inline">c_{i\mbox{ }\mbox{ }({\mu }_k)}</math> , must be answered.  For this we revisit the cases for boundary conditions.
+. ''Dirichlet boundary condition: ''Before imposing the boundary conditions at  <math display="inline">x=0</math> or  <math display="inline">y=0</math>, the nodal values along the two axes should be taken out from the vectors in null space of  [[Image:draft_Samper_725107104-image261.png|18px]] (Figure 4-b).  This may be accomplished by defining a suitable transformation matrix as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k=</math><math>\boldsymbol{T}_0^k{\phi }_{\mbox{ }\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (66)
+|}
+Here  <math display="inline">\boldsymbol{T}_0^k</math> contains the proportionality factors,  <math display="inline">{\mu }_1</math> and  <math display="inline">{\mu }_2</math> , and selects the nodal values of  [[Image:draft_Samper_725107104-image266.png|30px]] at boundaries.  Here again the dimension of  [[Image:draft_Samper_725107104-image267.png|6px]] depends on the number of the boundary patches considered along two axes.  Only few patches are practically needed for this purpose.  The prescribed boundary values in this case are as (32) instead of which, for simplicity, we use  <math display="inline">{\boldsymbol{\overline{u}}}_{prs}</math>.  Therefore
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{prs}=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left(\sum c_{i\mbox{ }\mbox{ }({\mu }_k)}^k\boldsymbol{v}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k\right)d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (67)
+|}
+We assume that the coefficient value associated with a value for <math display="inline">{\mu }_k</math> is proportional to the projection of  [[Image:draft_Samper_725107104-image271.png|24px]] on  [[Image:draft_Samper_725107104-image272.png|30px]] .  Such a proportionality assumption, which is consistent with conventional mathematical transformations, is applied in most general form as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| [[Image:draft_Samper_725107104-image273.png|126px]]
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (68)
+|}
+Where  [[Image:draft_Samper_725107104-image274.png|12px]] is a matrix to be determined so that Equation (68) holds for all values of  <math display="inline">{\mu }_k</math> .  To evaluate  [[Image:draft_Samper_725107104-image274.png|12px]] we substitute (68) in (67) to obtain
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{prs}=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k{\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right]\boldsymbol{R}{\boldsymbol{\overline{u}}}_{prs}d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (69)
+|}
+Since  [[Image:draft_Samper_725107104-image274.png|12px]] is assumed to be independent of  <math display="inline">{\mu }_k</math> it can be taken out of the integral
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{prs}=\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k{\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right]d{\mu }_k\right)\boldsymbol{R}\mbox{ }{\boldsymbol{\overline{u}}}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (70)
+|}
+Above relation implies that
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{R}={\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k{\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right]d{\mu }_k\right)}^{-1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (71)
+|}
+Our experience shows that  <math display="inline">\boldsymbol{R}</math> is well-conditioned and its inverse exists as long as enough integration points are used.  The number of integration points required is dependent on the number of prescribed values, i.e. the length of  <math display="inline">{\boldsymbol{\overline{u}}}_{prs}</math> .  The reader will note that vectors  <math display="inline">\boldsymbol{v}_i</math> are not orthogonal and thus increasing the number of integration points helps to obtain a full rank matrix.
+Having found the proportionality matrix  [[Image:draft_Samper_725107104-image274.png|12px]] , the finite element solution corresponding to the correction due to the presence of boundary is evaluated from (65) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \left\{{\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\boldsymbol{R}{\boldsymbol{\overline{u}}}_{prs}\right\}\boldsymbol{w}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k\right]\mbox{ }d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (72)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p=\left[\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left\{\sum \boldsymbol{w}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k{\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right\}d{\mu }_k\right]\boldsymbol{R}\mbox{ }{\boldsymbol{\overline{u}}}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (73)
+|}
+The rest of the procedure is construction of general asymptotic finite element solution as mentioned before.
+.  ''Neumann boundary condition: ''For this type of boundary conditions we revisit relation (37) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\left(\int_{\Omega^{'}_{bnd}} \boldsymbol{B}_{bnd}^T \boldsymbol{D} [\boldsymbol{B}_{bnd} \boldsymbol{B}_{1} ] d\Omega  \right) {\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd} =\boldsymbol{F}_{prs} </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (74)
+|}
+where  <math display="inline">\boldsymbol{B}_{bnd}</math> denotes a partition of  <math display="inline">\boldsymbol{B}</math> corresponding to the set of nodes at boundary on which  <math display="inline">\boldsymbol{F}_{prs}</math> is applied and  <math display="inline">\boldsymbol{B}_1</math> is the partition pertinent to nodes the shape function of which have projection to the nodes of  <math display="inline">\boldsymbol{B}_{bnd}</math>. Consequently  <math display="inline">{\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd}</math> in (74) represents a partition of  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> corresponding to nodes at or near boundary.  More briefly
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{K}_{bnd} {\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd} =  \boldsymbol{F}_{prs}  </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (75)
+|}
+Here  <math display="inline">\boldsymbol{K}_{bnd}</math> is a rectangular matrix of coefficients. We note that since in (74) just elements with at least one node on boundaries are contributing, a suitable transformation matrix as below should be defined
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}={\boldsymbol{T^{{'}}}}_0^{\mbox{ }\mbox{ }k}\boldsymbol{\phi }_{\mbox{ }\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (76)
+|}
+Where  <math display="inline">{\boldsymbol{T^{{'}}}}_0^{\mbox{ }\mbox{ }k}</math> contains the proportionality factors  <math display="inline">{\mu }_1</math> and  <math display="inline">{\mu }_2</math> and selects the nodal values of  <math>\boldsymbol{\phi }_{\mbox{ }\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math> at or near boundaries.
+Noting that a suitable subset of  <math display="inline">{\boldsymbol{w^{{k}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math> as  <math display="inline">{\boldsymbol{v^{{'k}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math>, similar to  <math display="inline">{\left({\boldsymbol{\overline{u}}}_{h0}^p\right)}_{bnd}</math> which is a subset of  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> , must be used then substitution of (65) in (75) leads to
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{K}_{bnd}\sum_{k=1}^2 {\int }_{\mbox{ }0}^{\mbox{ }1} \left(
+\sum c^{k}_{i(\mu_k)}\boldsymbol{v}^{\prime k}_{i({\mu }_k)}
+\right) d\mu_k= \boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (77)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\sum_{k=1}^2 {\int }_{\mbox{ }0}^{\mbox{ }1}\left(
+\sum c^{k}_{i(\mu_k)}\boldsymbol{K}_{bnd}\boldsymbol{v}^{\prime k}_{i({\mu }_k)}
+\right) d\mu_k= \boldsymbol{F}_{prs}   </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (78)
+|}
+In above   <math display="inline">\boldsymbol{K}_{bnd}\boldsymbol{v}^{\prime k}_{i({\mu }_k)}</math> represents the contribution of  <math display="inline">{\boldsymbol{w^{{k}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math> to nodal boundary resultants of tractions.  Therefore, similar to the first case, we assume that the coefficient value associated with a value for <math display="inline">\mu_k</math> is proportional to the projection of  <math display="inline">\boldsymbol{K}_{bnd}{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}</math> on  <math display="inline">\boldsymbol{F}_{prs}   </math>.  Hence
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>c_{i\mbox{ }\mbox{ }\mbox{ }({\mu }_k)}^k={\left(\boldsymbol{K}_{bnd}{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\boldsymbol{R^{{'}}}\mbox{ }\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (79)
+|}
+Where, again,  <math display="inline">\boldsymbol{R^{{'}}}</math> is to be determined so that Equation (79) holds for all values of  <math display="inline">{\mu }_k</math> .  To evaluate  <math display="inline">\boldsymbol{R^{{'}}}</math> we substitute (79) in (78) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{K}_{bnd}{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}{\left(\boldsymbol{K}_{bnd}{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]\boldsymbol{R^{{'}}}\mbox{ }\boldsymbol{F}_{prs}d{\mu }_k=</math><math>\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (80)
+|}
+and thus
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{K}_{bnd}\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum {\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}{\left({\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]d{\mu }_k\right)\boldsymbol{K}_{bnd}^T\boldsymbol{R^{{'}}}\mbox{ }\boldsymbol{F}_{prs}=</math><math>\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (81)
+|}
+Which implies that
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\textbf{R}'=\left[\textbf{K}_{bnd} \left(\sum _{k=1}^{2}\int_{0}^{1}\left[\sum \textbf{v}_{i(\mu _{k})}^{'k} (\textbf{v}_{i(\mu _{k})}^{'k})^{T}  \right]d\mu _{k}   \right)\textbf{K}_{bnd}^{T} \right]^{-1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (82)
+|}
+The finite element solution corresponding to the correction due to the presence of boundary is then evaluated from (65) as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math> \bar{\bf u}_{h0}^{p} =\sum _{k=1}^{2}\int _{\, 0}^{\, 1}\left(\sum \textbf{w}_{\, i \, (\mu _{k} )}^{k} \left[\left(\textbf{K}_{bnd} \textbf{v}_{i \, (\mu _{k} )}^{'k} \right)^{T} \textbf{R}'\; \textbf{F}_{prs} \right] \right)d\mu _{k} </math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (83)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p=\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{w}_{\mbox{ }i\mbox{ }\mbox{ }({\mu }_k)}^k{\left({\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]d{\mu }_k\right)\boldsymbol{K}_{bnd}^T\boldsymbol{R^{{'}}}\mbox{ }\boldsymbol{F}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (84)
+|}
+''3.  General mixed boundary condition.  ''To obtain a general relation between nodal resultants of tractions and associative nodal values, one may use (79) in (67).  Therefore by dropping the subscripts for  <math display="inline">\boldsymbol{F}</math> and [[Image:draft_Samper_725107104-image217.png|6px]] we have
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{\overline{u}}=\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{i\mbox{ }\mbox{ }\mbox{ }({\mu }_k)}^k{\left(\boldsymbol{K}_{bnd}{\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]\boldsymbol{R^{{'}}}\mbox{ }\boldsymbol{F}\mbox{ }d{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (85)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math> \bar{\bf u}=\sum _{k=1}^{2}\left(\int _{\, 0}^{\, 1}\left[\sum \textbf{v}_{i \, \, (\mu _{k} )}^{k} \left(\textbf{v}_{i \, (\mu _{k} )}^{'k} \right)^{T}  \right]d\mu _{k}  \right) \, \textbf{K}_{bnd}^{T} \textbf{R}'\; \textbf{F}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (86)
+|}
+in flexibility form or in stiffness form as
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{F}={\left[\sum_{k=1}^2\left({\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{i\mbox{ }\mbox{ }\mbox{ }({\mu }_k)}^k{\left({\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]d{\mu }_k\right)\mbox{ }\boldsymbol{K}_{bnd}^T\boldsymbol{R^{{'}}}\right]}^{-1}\boldsymbol{\overline{u}}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (87)
+|}
+Equation (87) allows us to use any mixed boundary conditions similar to a standard finite element solution.  Of course the associative  <math display="inline">{\boldsymbol{\overline{u}}}_{h0}^p</math> will be as (84) while dropping the subscript of  <math display="inline">\boldsymbol{F}</math> .
+It may be noticed that we could also use (73) in (74) for relating  <math display="inline">\boldsymbol{F}</math> to  <math display="inline">\boldsymbol{\overline u}</math> in one set of equations.  In that case <math display="inline">\boldsymbol{w}</math> should be replaced by  <math display="inline">\boldsymbol{v^{{'}}}</math> and therefore;
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{K}_{bnd}\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[{\sum {\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right]d{\mu }_k\right)\boldsymbol{R}\mbox{ }\boldsymbol{\overline{u}}=</math><math>\boldsymbol{F}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (88)
+|}
+It follows that theoretically the following relation must exist between  [[Image:draft_Samper_725107104-image274.png|12px]] and
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
 {|
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1157.svg|71px]]
+| <math>\boldsymbol{R^{{'}}}</math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1158.svg|center|80px]]
+| <math>\boldsymbol{K}_{bnd}\left(\sum_{k=1}^2{\int }_{\mbox{ }0}^{\mbox{ }1}\left[{\sum {\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\left(\boldsymbol{v}_{i\mbox{ }\mbox{ }({\mu }_k)}^k\right)}^T\right]d{\mu }_k\right)\boldsymbol{R}=</math><math>{\left[\sum_{k=1}^2\left({\int }_{\mbox{ }0}^{\mbox{ }1}\left[\sum \boldsymbol{v}_{i\mbox{ }\mbox{ }\mbox{ }({\mu }_k)}^k{\left({\boldsymbol{v^{{'}}}}_{i\mbox{ }\mbox{ }({\mu }_k)}^{\mbox{ }k}\right)}^T\right]d{\mu }_k\right)\boldsymbol{K}_{bnd}^T\boldsymbol{R^{{'}}}\right]}^{-1}</math>
 |}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (89)
+|}
+Of course from a numerical point of view there may be some discrepancy between the two relations.
+It is noteworthy to mention that the boundary conditions specified should be consistent with assumption made for the decay property, i.e.  <math display="inline">{\boldsymbol{\overline{u}}}_{prs}\rightarrow 0</math> or  <math display="inline">\boldsymbol{F}_{prs}\rightarrow 0</math> when  <math display="inline">\frac{d_x}{L_x}\rightarrow \infty </math> and   <math display="inline">\frac{d_y}{L_y}\rightarrow \infty </math> which of course is the case for our studies.
+If asymptotic finite element solution near a corner of a domain is sought, the decay condition  <math display="inline">\frac{d_z}{L_z}\rightarrow \infty </math> is also added to the above conditions.  To satisfy such an additional decay condition a corresponding proportionality property, with <math display="inline">{\mu }_3</math> as a parameter, between nodal values along  <math display="inline">z</math> direction must be considered.  The procedure leads to finding a characteristic equation with three parameters
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\vert {\boldsymbol{Q^{{'}}}}_{\left({\mu }_\boldsymbol{1},{\mu }_2,{\mu }_3\right)}\vert =</math><math>f({\mu }_1,{\mu }_2,{\mu }_3)=0</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (90)
+|}
+form which one parameter is determined by assuming values for the other two.  We define
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| {|
+|-
+| <math>{\mu }_3=h({\mu }_1,{\mu }_2)</math>
+| <math>h=\left\{g\mbox{ }\mbox{ }\mbox{ }\vert \mbox{ }\mbox{ }f({\mu }_1,{\mu }_2,g)=\right. </math><math>\left. 0,\mbox{ }\mbox{ }\mbox{ }\vert g\vert \leq 1,\mbox{ }\mbox{ }\vert {\mu }_1\vert \leq 1,\mbox{ }\mbox{ }\vert {\mu }_2\vert \leq 1\right\}</math>
+|}
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (91)
+|}
+Similar expressions may be written by replacing the parameters in turn.  Having found the null space of  <math display="inline">\boldsymbol{Q^{{'}}}</math> for each set, the independent nodal values of the patch is evaluated as
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math>{\boldsymbol{\overline{u}}}_I={\int }_{\mbox{ }0}^{\mbox{ }1}{\int }_{\mbox{ }0}^{\mbox{ }1}\sum c_i{}_{\mbox{ }({\mu }_1,{\mu }_2)}{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_1,{\mu }_2)}d{\mu }_1d{\mu }_2+</math><math>{\int }_{\mbox{ }0}^{\mbox{ }1}{\int }_{\mbox{ }0}^{\mbox{ }1}\sum c_i{}_{\mbox{ }({\mu }_1,{\mu }_3)}{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_1,{\mu }_3)}d{\mu }_1d{\mu }_3+</math><math>{\int }_{\mbox{ }0}^{\mbox{ }1}{\int }_{\mbox{ }0}^{\mbox{ }1}\sum c_i{}_{\mbox{ }({\mu }_2,{\mu }_3)}{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_2,{\mu }_3)}d{\mu }_2d{\mu }_3</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (92)
+|}
+or
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_I=\sum_{k=1}^2\sum_{n=k+1}^3{\int }_{\mbox{ }0}^{\mbox{ }1}{\int }_{\mbox{ }0}^{\mbox{ }1}\sum c_i{}_{\mbox{ }({\mu }_k,{\mu }_n)}{\phi }_{\mbox{ }\mbox{ }i\mbox{ }({\mu }_k,{\mu }_n)}d{\mu }_nd{\mu }_k</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (93)
+|}
+We have omitted the superscripts for convenience. Finally the formulation leads to
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>{\boldsymbol{\overline{u}}}_{h0}^p=\left[\sum_{k=1}^2\sum_{n=k+1}^3{\int }_{\mbox{ }0}^{\mbox{ }1}{\int }_{\mbox{ }0}^{\mbox{ }1}\sum \boldsymbol{w}_{i\mbox{ }({\mu }_k,{\mu }_n)}\boldsymbol{v}_{i\mbox{ }({\mu }_k,{\mu }_n)}^Td{\mu }_nd{\mu }_k\right]\boldsymbol{R}\mbox{ }{\boldsymbol{\overline{u}}}_{prs}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (94)
+|}
+Where
+{| class="formulaSCP" style="width: 100%; text-align: center;"
+|-
+|
+{| style="text-align: center; margin:auto;"
+|-
+| <math>\boldsymbol{R}={\left[\sum_{k=1}^2\sum_{n=k+1}^3{\int }_0^1{\int }_0^1\sum \boldsymbol{v}_{i\mbox{ }({\mu }_k,{\mu }_n)}\boldsymbol{v}_{i\mbox{ }({\mu }_k,{\mu }_n)}^Td{\mu }_nd{\mu }_k\right]}^{-1}</math>
+|}
+| style="width: 5px;text-align: right;white-space: nowrap;" | (95)
+|}
+Other relations may be found by analogy.
+It is observed that a finite element solution is achievable for Equation (49) with boundary conditions specified in (50) and (51).  The accuracy obtained depends upon the number of boundary patches nominated for satisfaction of boundary conditions at  <math display="inline">x=0</math> ,  <math display="inline">y=0</math> and  <math display="inline">z=0</math> .  In the limit, when the number of boundary patches tends to infinity, the solution must approach to a unique solution since the integral equation (49) and all boundary conditions are satisfied in a finite element sense. Our experience shows that few boundary patches are needed to obtain a reasonable converged solution. In a following sample problem we examine the convergence of the solution for a two dimensional case.
+=4. NUMERICAL PRESENTATION=
+''A sample problem ''
+Apart from assessment of error estimators, the method described above may be used to find a solution to a homogeneous differential equation like following Poisson’s equation
 {|
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1159.svg|71px]]
+| <math>{\nabla }^2u=0</math>  <math>x\geq 0</math> and    <math>y\geq 0</math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1160.svg|center|104px]]
 |}
+with exact solution, for instance, as
+<math display="inline">u_{ex}=e^{-(x+y)}cos(x-y)</math>
+It may be seen that the decay condition exists in the exact solution when  <math display="inline">x\rightarrow \infty </math> or  <math display="inline">y\rightarrow \infty </math> .  We try to give a finite element solution to such a problem using the proposed approach.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="text-align: center;vertical-align: top;width: 34%;"| [[Image:draft_Samper_725107104-image437.png|168px]]
+|  style="vertical-align: top;width: 34%;"| [[Image:draft_Samper_725107104-image438.png|168px]]
+|  style="text-align: center;vertical-align: top;width: 34%;"|             [[Image:draft_Samper_725107104-image439.png|120px]]
+|-
+|  style="text-align: center;vertical-align: top;"|     (a)
+|  style="text-align: center;vertical-align: top;"|     (b)
+|  style="text-align: center;vertical-align: top;"|(c)
+|}
+Figure 5.  Patch used for sample problem; (a) Union Jack pattern, (b) Four square elements as the basic patch, (c) A square element as the basic patch.
+''Solution with triangular elements''
+A mesh of infinite dimensions is considered with a repeatable pattern in which the smallest unit is of form as the patch shown in Figure 5-a.   After assembling the stiffness matrix of a macro patch of  <math display="inline">2\times 2</math> from the smallest unit, condensing the matrix for independent nodal values and also imposing the proportionality conditions similar to relations (52), the following expression for  <math display="inline">\boldsymbol{Q^{{'}}}</math> is obtained
+<math>\boldsymbol{Q^{{'}}}\left({\mu }_1,{\mu }_2\right)=</math><math>\left[\begin{array}{ccc}
+& -1-\frac{1}{{\mu }_1} & -1-\frac{1}{{\mu }_2}\\
+-1-{\mu }_1 & \frac{7}{2}-\frac{1}{4{\mu }_2}-\frac{{\mu }_2}{4} & -\frac{1}{4}-\frac{{\mu }_1}{4}-\frac{1}{4{\mu }_2}-\frac{{\mu }_1}{4{\mu }_2}\\
+-1-{\mu }_2 & -\frac{1}{4}-\frac{1}{4{\mu }_1}-\frac{{\mu }_2}{4}-\frac{{\mu }_2}{4{\mu }_1} & \frac{7}{2}-\frac{1}{4{\mu }_1}-\frac{{\mu }_1}{4}
+\end{array}\right]</math>
+The determinant of  <math display="inline">\boldsymbol{Q^{{'}}}</math> is set to be zero and thus the characteristic equation is obtained as
+<math>f\left({\mu }_1,{\mu }_2\right)=33+\frac{1}{4}\left({\mu }_1^2+\right. </math><math>\left. {\mu }_2^2+{\mu }_1^{-2}+{\mu }_2^{-2}\right)-</math><math>8\left({\mu }_1+{\mu }_2+{\mu }_1^{-1}+{\mu }_2^{-1}\right)-</math><math>\frac{1}{2}\left({\mu }_1{\mu }_2+{\mu }_1^{-1}{\mu }_2+\right. </math><math>\left. {\mu }_1{\mu }_2^{-1}+{\mu }_1^{-1}{\mu }_2^{-1}\right)=</math><math>0</math>
+We evaluate <math display="inline">{\mu }_2</math> in terms of  <math display="inline">{\mu }_1</math>
 {|
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1161.svg|104px]]
+| <math>{\mu }_2=g_1({\mu }_1)=\frac{1}{2}\left(16+{\mu }_1^{-1}- 8{\mu }_1^{-\frac{1}{2}}-8{\mu }_1^{\frac{1}{2}}+ {\mu }_1-{\mu }_1^{-1}\left(-1+5{\mu }_1^{\frac{1}{2}}- 5{\mu }_1+{\mu }_1^{\frac{3}{2}}\right){\left(1-6{\mu }_1^{\frac{1}{2}}+{\mu }_1\right)}^{\frac{1}{2}}\right)</math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1162.svg|center|113px]]
 |}
+<math>{\mu }_2=g_2\left({\mu }_1\right)=\frac{1}{2}\left(16+ {\mu }_1^{-1}+8{\mu }_1^{-\frac{1}{2}}+8{\mu }_1^{\frac{1}{2}}+{\mu }_1-{\mu }_1^{-1}\left(1+5{\mu }_1^{\frac{1}{2}}+ 5{\mu }_1+{\mu }_1^{\frac{3}{2}}\right){\left(1+6{\mu }_1^{\frac{1}{2}}+{\mu }_1\right)}^{\frac{1}{2}}\right)</math>
+<math>{\mu }_2=g_3\left({\mu }_1\right)=g_1{\left({\mu }_1\right)}^{-1},\mbox{ }{\mu }_2=g_4\left({\mu }_2\right)=g_2{\left({\mu }_1\right)}^{-1}</math>
+Similar expressions may be derived for  <math display="inline">{\mu }_1</math> in terms of  <math display="inline">{\mu }_2</math> . Generally, such explicit relations are not necessary and what needed is just choosing a value for  <math display="inline">{\mu }_1</math> , as an integration point in the feasible domain, and solving an algebraic equation for its roots. Nevertheless, to give some insight to the problem the explicit relations are given here. Now we define
+{|
+|-
+| <math>{\left({\mu }_2\right)}_1=h_1({\mu }_1)</math>
+| <math>h_1=\left\{g\mbox{ }\vert \mbox{ }\mbox{ }g\in \left\{g_1,g_3\right\}\mbox{ },\mbox{ }\vert g\vert \leq 1\right\}</math>
+|}
+,
 {|
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1163.svg|73px]]
+| <math>{\left({\mu }_2\right)}_2=h_2({\mu }_1)</math>
-| [[Image:draft_Samper_725107104-picture-x0000_i1164.svg|center|72px]]
+| <math>h_2=\left\{g\mbox{ }\vert \mbox{ }\mbox{ }g\in \left\{g_2,g_4\right\}\mbox{ },\mbox{ }\vert g\vert \leq 1\right\}</math>
 |}
-[[Image:draft_Samper_725107104-picture-x0000_i1165.svg|104px]]
-In order to study the boundary layer solution, variation of [[Image:draft_Samper_725107104-picture-x0000_i1166.svg|23px]] near corner is illustrated in Figure 2.  In this solution 10 patches are used along each axis.
-It may be noticed that it is not easy to give a direct finite element solution for comparison since no clear information is available for boundary conditions at the other two edges of the domain.  However, one may take either of two conditions given in (28-[[Image:draft_Samper_725107104-picture-x0000_i1167.svg|19px]] ) in an approximate manner, i.e. [[Image:draft_Samper_725107104-picture-x0000_i1168.svg|48px]] or [[Image:draft_Samper_725107104-picture-x0000_i1169.svg|59px]] .  This means that approximately either a homogeneous Dirichlet or Neumann boundary condition might be chosen for far ends of the domain.  Depending on the main boundary conditions near the corner, each case may have some drawbacks such as presence of singular points and etc.   However, if the domain, nominated for direct finite element solution, is chosen very large then the effects of such imperfections at far ends of the domain become negligible at the patches near the corner.   Our experience shows that use of Neumann boundary condition gives good results when the dimensions of the domain are very large.  Figure 3 depicts the variation of [[Image:draft_Samper_725107104-picture-x0000_i1170.svg|13px]] obtained form such approximate and direct finite element solution.  Comparison of the results with those shown in Figure 2, for the proposed method, illustrates the similarities between the two sets (the difference between the contour presentations in the two figures is basically due to the difference between the levels of solutions after the first layers of the elements).
+when    <math display="inline">\vert {\mu }_1\vert \leq 1</math> .  Similar functions are defined, with analogy, when  <math display="inline">{\mu }_1</math> is evaluated in terms of  <math display="inline">{\mu }_2</math> .
-''Sample problem 2''
+We should note that generally all quantities may take complex values, however, as pointed out before we take real values for one parameter and calculate the other one which may be in the form of a complex pair.
-In this example we shall use regular pattern of triangles as Figure 4-a. Once again, the problem is to find [[Image:draft_Samper_725107104-picture-x0000_i1171.svg|23px]] near corner of the domain.  For the asymptotic finite element solution, the boundary condition to be imposed is of a Neumann form extracted from monomial[[Image:draft_Samper_725107104-picture-x0000_i1172.svg|27px]] as an exact solution.    For convenience in giving numeric details, we use a patch of two elements as the basic repeating unit of the mesh.
+To give a finite element solution, we take exact fluxes normal to boundaries at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> as tractions for the infinite domain to be solved.  We note that although relations (87) or (88) may be employed to first find the unknown nodal boundary values and then extend the result to the interior parts through Equation (73) , one may directly use Equation (84) instead.  Figure 6 depicts the variation of the finite element solution for  <math display="inline">u</math> obtained from the proposed method using 10 patches along each axis,  <math display="inline">x=0</math> and  <math display="inline">y=0</math> .  In order to show that the solution is close to an approximate finite element solution, in Figure 7 we give the results of a finite element solution on a domain of  <math display="inline">10\times 10</math> patches using exact fluxes for the other two boundaries.  Although we do not expect that the two sets of results be identical, due to the fact that the tractions used for the latter finite element solution do not represent the tractions of remaining domain in FEM sense, the results are very close (the reader may also note that since all the boundary conditions for the latter solution are in Neumann form one node of the domain should be restrained, the most remote one in this case, thus this constitutes another source for the difference between the two sets of results).   It can clearly be seen that the method proposed is able to give a finite element solution for a domain with unbounded dimension exhibiting a decay property.
-According to the node numbering in Figure 4-b, the nodal values of periodic finite element solution, evaluated from Equation (23-[[Image:draft_Samper_725107104-picture-x0000_i1173.svg|19px]] ) and conditions (24-[[Image:draft_Samper_725107104-picture-x0000_i1174.svg|19px]] ),  are as
-[[Image:draft_Samper_725107104-picture-x0000_i1175.svg|216px]]
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image440.png|288px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image441.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-The tractions are calculated from (33-[[Image:draft_Samper_725107104-picture-x0000_i1176.svg|19px]] ).  Then the prescribed nodal forces for three patches (here elements) along [[Image:draft_Samper_725107104-picture-x0000_i1177.svg|37px]] and[[Image:draft_Samper_725107104-picture-x0000_i1178.svg|37px]] are obtained as (see numbering of the nodes in Figure 4-b)
+Figure 6. Finite solution obtained from the proposed method for the sample problem using  <math display="inline">10\times 10</math> patches with '''Union Jack''' pattern in Fig. 5-a and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
-[[Image:draft_Samper_725107104-picture-x0000_i1179.svg|56px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1180.svg|137px]]
-Note that five nodes are nominated for satisfaction of the boundary conditions (nodes shown in dashed box in Figure 4-b) though twelve nodes are involved in the computation since seven nearby nodes are contributing to traction resultants at nominated nodes through the adjacent elements shown in shaded gray area in Figure 4-b.  Stiffness matrix of the patch is evaluated as
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 52%;"| [[Image:draft_Samper_725107104-image442.png|288px]]
+|  style="vertical-align: top;width: 49%;"| [[Image:draft_Samper_725107104-image443.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-{|
+Figure 7. Standard finite solution for the sample problem using  <math display="inline">10\times 10</math> patches with '''Union Jack''' pattern in Fig. 5-a and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> , and other two boundaries; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
+The accuracy of the solution may be increased or decreased by use of different number of patches along the two axes.  Figures 8 and 9 show the variation of  <math display="inline">u</math> over the domain using different number of patches along each direction.
+{| style="width: 100%;border-collapse: collapse;"
 |-
-| [[Image:draft_Samper_725107104-picture-x0000_i1181.svg|237px]]
+|  style="vertical-align: top;width: 52%;"| [[Image:draft_Samper_725107104-image444.png|288px]]
-| [[Image:draft_Samper_725107104-picture-x0000_i1182.svg|center|156px]]
+|  style="vertical-align: top;width: 49%;"| [[Image:draft_Samper_725107104-image445.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
 |}
-,      [[Image:draft_Samper_725107104-picture-x0000_i1183.svg|156px]]
-We apply the proportionality relations to obtain the following characteristic equation
-[[Image:draft_Samper_725107104-picture-x0000_i1184.svg|231px]]
+Figure 8. Finite solution obtained from the proposed method for the sample problem using  <math display="inline">8\times 8</math> patches with '''Union Jack''' pattern in Fig. 5-a and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> . The results are shown on  <math display="inline">10\times 10</math> patches; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
-Note that in this simple case no null space solution is required since just one independent DOF is remaining.  Now for an integration point as [[Image:draft_Samper_725107104-picture-x0000_i1185.svg|55px]] one may calculate [[Image:draft_Samper_725107104-picture-x0000_i1186.svg|147px]] form above characteristic equation.  We construct [[Image:draft_Samper_725107104-picture-x0000_i1187.svg|17px]] and [[Image:draft_Samper_725107104-picture-x0000_i1188.svg|19px]] for nodes 1 to 12, in Figure 4-b, from the proportionality factors as
-[[Image:draft_Samper_725107104-picture-x0000_i1189.svg|600px]]
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 52%;"| [[Image:draft_Samper_725107104-image446.png|276px]]
+|  style="vertical-align: top;width: 49%;"| [[Image:draft_Samper_725107104-image447.png|270px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-Also [[Image:draft_Samper_725107104-picture-x0000_i1190.svg|33px]] is evaluated, for DOFs shown in shaded area of Figure 6-b, as
+Figure 9. Finite solution obtained from the proposed method for the sample problem using  <math display="inline">6\times 6</math> patches with '''Union Jack''' pattern in Fig. 5-a and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> . The results are shown on  <math display="inline">10\times 10</math> patches; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
-[[Image:draft_Samper_725107104-picture-x0000_i1191.svg|600px]]
-Performing numerical integration in (82-[[Image:draft_Samper_725107104-picture-x0000_i1192.svg|19px]] ) for points selected for both [[Image:draft_Samper_725107104-picture-x0000_i1193.svg|19px]] and [[Image:draft_Samper_725107104-picture-x0000_i1194.svg|20px]] , matrix [[Image:draft_Samper_725107104-picture-x0000_i1195.svg|21px]] is evaluated as
+Figure 10 also show the convergence of the results in terms of the number of integration points and patches used.  In these figures the nodal values obtained for nodes 1 and 5 of the corner patch in Figure 4-b are plotted against the number of Gauss points and patches used.  At the first glimpse no uniform convergence is seen, however, more careful look into the results reveals that the two values for nodes 1 and 5 follow each other with constant difference.  The reason for such behavior lies in the rigid body motion of the whole sub-domain.  We notice that the boundary conditions used are of Neumann type and no explicit essential condition is used except for the decay condition the strength of which varies between different solutions.  The gradients and  relative values show excellent convergence for moderate number of Gauss points and patches as shown in Figure 10 for  <math display="inline">u_5-u_1</math> and  <math display="inline">{\sigma }_x</math> , the gradient at center of element 1-2-9 in Figure 5-a.
-[[Image:draft_Samper_725107104-picture-x0000_i1196.svg|600px]]
-Now we define[[Image:draft_Samper_725107104-picture-x0000_i1197.svg|20px]] for nodes 1, 2, 3, 5, 6 and 8 for example, through the following transformation
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"|
+[[File:Draft_Samper_725107104_8068_Fig10a.png]]
+|  style="vertical-align: top;width: 51%;"|
+[[File:Draft_Samper_725107104_4232_Fig10b.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-[[Image:draft_Samper_725107104-picture-x0000_i1198.svg|276px]] ,    [[Image:draft_Samper_725107104-picture-x0000_i1199.svg|173px]] ,   [[Image:draft_Samper_725107104-picture-x0000_i1200.svg|175px]]
-Boundary layer solution, [[Image:draft_Samper_725107104-picture-x0000_i1201.svg|23px]] , is obtained from Equation (84-[[Image:draft_Samper_725107104-picture-x0000_i1202.svg|19px]] )
+Figure 10.  Convergence of the solution with Neumann boundary condition-  <math display="inline">u_1</math> and  <math display="inline">u_5</math> are nodal values of node one and five in '''Union Jack''' Pattern of Fig 5-a for the corner patch of the domain, and  <math display="inline">\sigma </math> is the gradient at center of triangle 1-2-9; (a) Convergence with respect to the number of Gauss points in a  <math display="inline">10\times 10</math> macro patch, (b) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.
-[[Image:draft_Samper_725107104-picture-x0000_i1203.svg|125px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1204.svg|180px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1205.svg|180px]] , [[Image:draft_Samper_725107104-picture-x0000_i1206.svg|127px]]
-The asymptotic finite element solution for the patch is determined from [[Image:draft_Samper_725107104-picture-x0000_i1207.svg|144px]] where [[Image:draft_Samper_725107104-picture-x0000_i1208.svg|101px]] is calculated from (16-[[Image:draft_Samper_725107104-picture-x0000_i1209.svg|19px]] ).  For patch shown in Figure 4-b
-[[Image:draft_Samper_725107104-picture-x0000_i1210.svg|120px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1211.svg|169px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1212.svg|169px]] ,  [[Image:draft_Samper_725107104-picture-x0000_i1213.svg|112px]]
+We now examine the performance of the method when Dirichlet type of boundary conditions is of concern. The exact values are taken as prescribed boundary conditions along  <math display="inline">x=0</math> and  <math display="inline">y=0</math> . Equation (73) is directly employed to obtain the solution.  To give a finite element solution for comparison, here again a square domain is solved with the prescribed values, at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> , and Neumann boundary conditions at the other two edges.  The results of the two solutions are presented in Figures 11 and 12.  It can again be seen that the proposed method is able to produce the finite element solution.  Figure 13 depicts the convergence of the proposed solution, in terms of the number of integration points and patches used.  Unlike the case for Neumann boundary conditions, in this case uniform and fast convergence is seen in the results.
-We note that the values obtained for [[Image:draft_Samper_725107104-picture-x0000_i1214.svg|23px]] may still have slight rigid body motion since the decay condition has not been satisfied effectively and thus the values of [[Image:draft_Samper_725107104-picture-x0000_i1215.svg|17px]] are of high error except for the gradients at element level (see also the discussion for the sample problem in [1]).  The rest of the procedure for finding upper and lower bounds of effectivity index is straight forward and may be found in literature.
-It is worthy to make some insight for distribution of boundary layer solution.   Figure 5 shows the solution when twenty patches are used along each axis.  For comparison, we have included the results for [[Image:draft_Samper_725107104-picture-x0000_i1216.svg|45px]] larger patches constructed from [[Image:draft_Samper_725107104-picture-x0000_i1217.svg|33px]] patches with similar pattern.  It can be seen that the solution decays fast towards the interior parts. In another point of view, the solution suggests the use of homogenous Neumann boundary condition if an approximate direct FEM solution is to be applied on a finite domain.   When Neumann conditions are of concern at main boundaries, as is the case of this example, use of Neumann conditions at the other edges of the domain makes some difficulties for obtaining the solution.   However, to solve the problem one may restrain the most remote DOFs.  Of course, as discussed in the first example, this may introduce a sort of singularity in the domain.  Figure 6 show a standard finite element solution for the problem when Neumann boundary conditions have been used for far edges and the most remote DOF is restrained.  The results show reasonable agreement between the two solutions near the corner.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image448.png|288px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image447.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-==3. The error estimators==
-In this study we shall use some recovery based error estimators proven to be robust in two dimensional problems.  Apparently, robustness of a recovery based estimate lies in the quality of the recovered fields used and this means that the recovery should show superconvergent or semi-superconvergent property.  A recovery based error estimate is said to be “asymptotically exact” when the recovery procedure used is at least “asymptotically superconvergent”. We shall examine the performance of the recoveries when applied to three dimensional cases.
+Figure 11. Finite solution obtained from the proposed method for the sample problem using  <math display="inline">10\times 10</math> patches with '''Union Jack''' pattern in Fig. 5-a and Dirichlet boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
-The use of recovery based error estimates in robustness test is consistent with assumption made, in the first part, for smoothness of the exact solution.  Generally the test results may not be used for predicting the performance of the recovery method near singularities. Recent studies by Babuška ''et al'' on use of recovery based error estimates in problems with singularities show that pollution error is dominant over the whole domain [12-14] which may not be captured by use of local recovered information in the error estimation.  However, as we discussed in the introduction this deficiency is practically overcome in an adaptive procedure when the mesh becomes finer at singularities. Excessive refinement around singular points leads to fast reduction of the pollution error and thus even at points not very far from the singularities the assumptions made for smoothness of the solution will be valid.  Therefore if the estimated error shows reasonable asymptotic behavior the results of estimated error becomes more accurate while the mesh becomes finer.'' ''
-===3.1 Recoveries used===
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 52%;"| [[Image:draft_Samper_725107104-image449.png|294px]]
+|  style="vertical-align: top;width: 49%;"| [[Image:draft_Samper_725107104-image450.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-''Superconvergent Patch Recovery (SPR)''
-One of the most robust recovery methods was proposed by Zienkiewicz and Zhu in 1992.  Although at the time of its proposal some benchmark problems were used to demonstrate its performance, later comprehensive studies by Babuška and co-workers [9-11] proved its superiority over the other methods in terms of robustness and asymptotic convergence rate in two dimensional problems. It is noteworthy to mention that the studies in [9-11] show that none of the augmented forms of SPR, as some researches suggest [6-8] can compete with the original formulation.  In this paper we shall perform a comprehensive study on the performance of SPR when applied to three dimensional problems.
+Figure 12. Standard finite solution for the sample problem using  <math display="inline">10\times 10</math> patches with '''Union Jack''' pattern in Fig. 5-a and Dirichlet boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> , and Neumann condition at the other two boundaries; (a) Contour line presentation, (b) Three dimensional plot presentation with horizontal axes scaled with element size.
-The procedure of SPR may be found in many texts and papers [4, 5, 21-24].  The procedure starts with assuming a continuous gradient field, in the form of a polynomial, over a patch of elements built on a vertex node
-[[Image:draft_Samper_725107104-picture-x0000_i1218.svg|61px]] (2)
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"|
+[[File:Draft_Samper_725107104_4964_Fig13a.png]]
+|  style="vertical-align: top;width: 51%;"|
+[[File:Draft_Samper_725107104_2540_Fig13b.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-In which [[Image:draft_Samper_725107104-picture-x0000_i1219.svg|13px]] contains the monomials, in this paper a complete polynomial in a three dimensional space, and [[Image:draft_Samper_725107104-picture-x0000_i1220.svg|13px]] is a matrix of the unknown coefficients.  Minimization of a functional, as an error norm built from summation of errors between the enhanced field and the finite element solution at some sampling points, as
-[[Image:draft_Samper_725107104-picture-x0000_i1221.svg|129px]] (3)
+Figure 13.  Convergence of the solution with Dirichlet boundary condition-  <math display="inline">u_1</math> and  <math display="inline">u_5</math> are nodal values of node one and five in '''Union Jack''' Pattern of  Fig 5-a for the corner patch of the domain, and  <math display="inline">\sigma </math> is the gradient at center of triangle 1-2-9; (a) Convergence with respect to the number of Gauss points in a  <math display="inline">10\times 10</math> macro patch, (b) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.
-with respect to polynomial coefficients [[Image:draft_Samper_725107104-picture-x0000_i1222.svg|13px]] results to the following expression for local recovered stress/gradient field
-[[Image:draft_Samper_725107104-picture-x0000_i1223.svg|188px]] (4)
+''Solution with four node elements''
-The enhanced recovered gradient field for the whole domain is constructed by calculation of the above local field at patch point and interpolation of the nodal values obtained for the whole domain via shape functions used for the main unknowns [[Image:draft_Samper_725107104-picture-x0000_i1224.svg|19px]] .
+A mesh consisting of repeated patterns of a four square element as a macro-patch shown in Figure 5-b is considered for solution.   Applying the proportionality relations, as (52), the following expression for  the characteristic equation is obtained
-In this paper the sampling points used for three dimensional linear elements are shown in Figure 7.   For further details of the procedure the reader may consult to references [4, 5, 21-24].
+<math>f\left({\mu }_1,{\mu }_2\right)=\frac{799}{54}-\frac{259}{108}F_1+</math><math>\frac{25}{54}F_2-\frac{1}{108}F_3-\frac{49}{108}F_1^2+</math><math>\frac{1}{108}F_1F_2-\frac{1}{648}F_1F_3+\frac{1}{864}F_2^2+</math><math>\frac{1}{2592}F_1^4</math>
-''Recovery by Equilibrium in Patches (REP) ''
+Where
-Another superconvergent recovery method, REP,  was first proposed in [15] and then improved in [16]. In the latter research it was shown that although the original form of REP loses its robustness for patches with high aspect ratios, the improved version is as robust as SPR and can be considered as an alternative when the location of sampling points are not defined a priori.   The procedure is based on an integral form of equilibrium equation.  The idea is to find a local continuous field of gradient so that the equilibrium of patch is preserved, at least, approximately
+<math>F_1={\mu }_1+{\mu }_1^{-1}+{\mu }_2+{\mu }_2^{-1}</math> ,            <math>F_2={\mu }_1^2+{\mu }_1^{-2}+{\mu }_2^2+{\mu }_2^{-2}</math> ,        <math>F_3={\mu }_1^3+{\mu }_1^{-3}+{\mu }_2^3+{\mu }_2^{-3}</math>
-[[Image:draft_Samper_725107104-picture-x0000_i1225.svg|112px]] ,       [[Image:draft_Samper_725107104-picture-x0000_i1226.svg|116px]] (5)
+From above characteristic equation  <math display="inline">{\mu }_2</math> is to be calculated in terms of values given for  <math display="inline">{\mu }_1</math> and vice versa. Two forms of problem, i.e. with Neumann and Dirichlet conditions, are solved using the proposed method. Some results obtained from using different numbers of patches along the axes are given in Figures 14 and 15.  Once again, in order to compare the results with a standard finite element solution, a square domain is solved with the conditions mentioned previously, i.e. tractions from exact solution at all edges (see Figure 16).  Convergence of the solution for one node of the domain is illustrated in Figure 17 in terms of number of integration points and patches used. Similar discussion as given for Figure 10 is valid here and as is seen the relative quantities and the gradients exhibit excellent convergence.
-In above subscript [[Image:draft_Samper_725107104-picture-x0000_i1227.svg|16px]] denotes patch-wise view of notations, [[Image:draft_Samper_725107104-picture-x0000_i1228.svg|19px]] represent the matrix of stresses/gradients each of its components is expressed as (2) and [[Image:draft_Samper_725107104-picture-x0000_i1229.svg|20px]] is a vector of forces induced by finite element stress/gradient field maintaining the patch in equilibrium state.  The procedure resembles to the solution of a system, as small as a patch, using a set of piecewise continuous weight functions and a set of generalized coordinates as the polynomial terms.  The formulation can be written for strain or even displacement field.
-The original form of REP, proposed in [15], directly employs Equation (5) to construct enhanced field of gradients in a coupled form, while the improved version uses spectral form of (5) to construct the recovered gradient field for each gradient component separately.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image451.png|264px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image452.png|270px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-''Original form''
-In this version a least square procedure is used for reducing the difference between forces induced by finite element gradients and those induced by the recovered fields.  Hence the following functional
+Figure 14. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of '''four squares''' in Fig. 5-b and Dirichlet boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Solution for  <math display="inline">8\times 8</math> patches shown on a domain of <math display="inline">10\times 10</math> patches, (b) Solution for  <math display="inline">10\times 10</math> patches.
-[[Image:draft_Samper_725107104-picture-x0000_i1230.svg|289px]] (6)
-is minimized with respect to coefficients in [[Image:draft_Samper_725107104-picture-x0000_i1231.svg|19px]] leading to a coupled system of equations.  The rest of the procedure is similar to that of SPR mentioned above.  The reader may also consult to the references 15,16 and 24 for further explanation.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image453.png|276px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image454.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-''Improved form''
-In reference 16 it has been shown that the formulation given for original REP loses robustness for meshes with high aspect ratio. The sensitivity to aspect ration is effectively rediuced  if we write both the finite element and recovered stress vectors, [[Image:draft_Samper_725107104-picture-x0000_i1232.svg|21px]] and [[Image:draft_Samper_725107104-picture-x0000_i1233.svg|20px]] , as:
+Figure 15. Contour plots the finite solutions obtained from the proposed method for the sample problem using pattern of '''four squares''' in Fig. 5-b and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Solution for  <math display="inline">8\times 8</math> patches shown on a domain of <math display="inline">10\times 10</math> patches, (b) Solution for  <math display="inline">10\times 10</math> patches.
-[[Image:draft_Samper_725107104-picture-x0000_i1234.svg|273px]] (7)
-and substitute in Equation (5) to obtain
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image455.png|276px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image449.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-[[Image:draft_Samper_725107104-picture-x0000_i1235.svg|221px]] (8)
-and split the equilibrium equations accordingly.  The procedure continues with minimization of a series of functionals as
+Figure 16. Contour plots for standard finite solutions for the sample problem using  <math display="inline">10\times 10</math> pattern of '''four squares''' in Fig. 5-b and boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Neumann condition, (b) Dirichlet condition.
-[[Image:draft_Samper_725107104-picture-x0000_i1236.svg|600px]] (9)
-with respect to coefficients in [[Image:draft_Samper_725107104-picture-x0000_i1237.svg|19px]] .   The rest of procedure is again similar to SPR as mentioned before.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_1191_Fig17a.png]]
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_7634_Fig17b.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|-
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_2837_Fig17c.png]]
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_1279_Fig17d.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(c)
+|  style="text-align: center;vertical-align: top;"|(d)
+|}
-''Enhanced gradients at boundaries''
-Because in this paper assessment of the recoveries at boundaries of the domain is aimed, it is necessary to specify the way that nodal values of recovered gradients are calculated.  The reader may note that there are several ways for gradient calculation at boundaries.  An efficient way which was suggested in [23, 24] is using the information from patches nearby the boundary node.  This means that one can directly use the smooth field obtained in (2) for patches near boundary and substitute the coordinates (or local normalized coordinates) of the boundary node in the polynomial terms.  Since a node at the boundary may generally fall in several patches, especially in three dimensional domains, it is preferred to obtain the information from patch of closest internal nodal point.    However, this makes some ambiguity when unstructured meshes are used since the distances between a boundary node and the nodes around it may not be so different.  Therefore, here, we employ an averaging approach.
+Figure 17.  Convergence of the solution with Neumann and Dirichlet boundary condition-  <math display="inline">u_1</math> and  <math display="inline">u_5</math> are nodal values of node one and five in Fig 5-b for the corner patch of the domain, and  <math display="inline">\sigma </math> is the gradient at center of element 1-2-9-8; (a)/(c) Convergence of solution with respect to the number of Gauss points in a  <math display="inline">10\times 10</math> macro patch, (b)/(d) Convergence of solution with respect to the number of patches used along each direction while the number of Gauss points is fixed to 30 for all solutions.
-==4.  TEST RSULTS==
-In this section we shall present the results of the tests on some patterns of linear elements for three dimensional heat conduction and elasticity problems.  The results for elasticity problems are taken for Poisson’s ratio equal to 0.3.
-The three dimensional patterns of tetrahedral elements are recognized from the two dimensional patterns appearing on the faces of the cuboid.  Four well known patterns of triangles in two dimensional problems as shown in Figure 8 are used for naming the three dimensional patterns (see Figure 9).  For example RRR, according to this terminology, represents a 3D pattern having Regular arrangements of elements at three faces.  Since there will be two cases which can be named as CCC, i.e. the cuboids with Criss-cross and Chevron faces, we recognize the former as III.
+It is interesting to note that for solution using square elements, the smallest repeatable pattern may be considered as one element.  Therefore, it is worthwhile to compare the results from such a solution with those obtained above.
-The results of the tests for tetrahedral elements are summarized in Tables 1 to 4 for different aspect ratios.  We have also included the results of linear brick elements in Table 5.
+If the proportionality conditions are written for one element, just one relation is obtained for  <math display="inline">\boldsymbol{Q^{{'}}}</math> which is identical to the characteristic equation
-All Tables contain the lower and upper bounds as well as the robustness values for SPR and the two versions of REP.   In each table the results of the test at interior parts with periodicity in all directions, at flat boundaries with periodicity in two directions, and at kinked boundaries with periodicity in one direction, are listed.  Except for the tests at interior parts, two sets of results pertaining to the type of boundary conditions, i.e. Neumann or Dirichlet type, are given in each case.
+<math>f\left({\mu }_1,{\mu }_2\right)=\frac{1}{3}\left[8-\right. </math><math>\left. \left({\mu }_1+{\mu }_1^{-1}+{\mu }_2+{\mu }_2^{-1}+\right. \right. </math><math>\left. \left. {\mu }_1{\mu }_2+{\mu }_1{\mu }_2^{-1}+\right. \right. </math><math>\left. \left. {\mu }_1^{-1}{\mu }_2+{\mu }_1^{-1}{\mu }_2^{-1}\right)\right]</math>
-The results for patches at corners, where the periodicity is lost in all directions, with Neumann type of condition are given in Table 6. We have not reported the results of patches at corners with Dirichlet condition since the frequency of happening such situations in practical engineering problems is much smaller than the other cases presented.
+Figures 18 to 20 show the results obtained.  Here  <math display="inline">20\times 20</math> patches are used in order that the results be comparable with those of Figures 14,15 and 17.  Study on the convergence of the solution, while comparing to the solution with  <math display="inline">2\times 2</math> elements, demonstrates a slower rate of convergence in terms of the number of Gauss points and a faster rate in terms of the number of patches used (Figure 20).
-When the patch is located at interior of the domain, the tables clearly show that all recovery methods give superconvergent results, for all patterns except for the CCC.  For high aspect ratios of CCC pattern, however, REP in its original form fails to give superconvergent results while the improved version of REP and SPR still give satisfactory results but of course not fully superconvergent.  Superconvergence is also seen for brick elements for patches at interior parts.
-Interesting to note is that when patches are located at boundaries, none of the recovery methods are able to produce full superconvergent results.  However, the results for upper and lower values are within satisfactory ranges.  Exception may be seen for III pattern, i.e. a cuboid with Criss-cross pattern at all faces, when used near Neumann boundaries for which the lower bound values are less than the other cases.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image456.png|264px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image457.png|270px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-Similar conclusions may be made from the results obtained for elasticity problems.  Tables  7 to 12 show the robustness obtained for this class of problems.  Although the robustness indices are still in an acceptable range, they are of higher values than those obtained for heat conduction problems.  There is also an exception for results of brick elements, Table 11, in which full superconvergent effect is still seen from application of SPR similar to the heat cases. Comparison of all lower bounds for tetrahedral elements shows, again, that those of III pattern are of less value.
-The main reason for loss of superconvergent effect lies in the definition of the recoveries near boundaries.  This has also been quoted in studies of reference [20] performed for two dimensional problems with flat boundaries.  In fact, the minimizations performed in each of the recoveries leads to best answers near center of the patch, due to concentration of information around the patch point.  This inherently means that, even for interior parts of the mesh, for some types of elements the polynomial fitted in the recovery does not thoroughly capture the exact field. Nevertheless, in such cases it usually happens that the polynomial passes the exact solution at the patch node or at vicinity of it and since the recovered field is constructed by interpolation of nodal values the final result may show full on semi superconvergent effect.
+Figure 18. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of '''one square''' element in Fig. 5-c and Dirichlet boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Solution for  <math display="inline">16\times 16</math> patches shown on a domain of <math display="inline">20\times 20</math> patches, (b) Solution for  <math display="inline">20\times 20</math> patches.
-Another reason for loss of superconvergent effect seems to be the presence of spurious boundary layer modes.  This kind of undesirable solution usually vanishes within a distance of few elements length.  The affected distance is obviously dependent on the type of elements used. This latter reason seems to be less pronounced than the first one since our experience shows that the amplitude of the spurious modes is small compared to the nodal values of asymptotic finite element solution. However, this observation confirms that producing enhanced gradients from local information at elements on the boundaries may not lead to a full superconvergent effect especially when the whole patch, used in the recovery procedure, falls inside the affected region.  In this respect, the behavior is similar to cases with presence of a singularity in the domain [12-14].  But we note that in the present cases the exact solutions are smooth.
-From above discussion it may be concluded that two remedies may be thought for loss of superconvergent effect.  One is to use information from inner layers, beyond the first layer, and extrapolate the effect to the boundaries while devising some means for preventing lack of fitness in the least square procedure.  This may be performed by extrapolating the recovered nodal values obtained for patch nodes at neighborhood of the boundary node.  Another remedy, which should be considered simultaneously with the first one, is capturing the boundary layer modes and removing them from the enhanced solution which is not an easy task but still possible by considering that in this paper we have introduced a suitable tool for study of behavior of boundary layers.  Both remedies are beyond of the scopes of this paper and we shall refer to these in subsequent contributions.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image458.png|276px]]
+|  style="vertical-align: top;width: 51%;"| [[Image:draft_Samper_725107104-image459.png|282px]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|}
-As a final conclusion from comparison between SPR and improved REP applied to heat and elasticity problems, it may be stated that although full superconvergent property is lost for tetrahedral elements, results from both recoveries exhibit similar robustness.  However for brick elements, both methods produce full superconvergent gradients in an asymptotic sense for heat problems, but for elasticity problems, just SPR is able to produce full superconvergent effect.
-''Complementary results ''
+Figure 19. Contour plots for the finite solutions obtained from the proposed method for the sample problem using pattern of '''one square''' element in Fig. 5-c and Neumann boundary condition at  <math display="inline">x=0</math> and  <math display="inline">y=0</math> ; (a) Solution for  <math display="inline">16\times 16</math> patches shown on a domain of <math display="inline">20\times 20</math> patches, (b) Solution for  <math display="inline">20\times 20</math> patches.
-As is seen from tests results for interior parts, full superconvergent effect is lost for some mesh configuration.  Here we shall investigate the possibility of losing the effect by producing mesh configuration form combination of patterns at three main faces of a cuboid. The aspect ratios considered are identical to those given in the first columns of Tables 1 to 12.  To give more insight to the problem we present the results of robustness indices with respect to average minimum angle between the faces of tetrahedral elements in a patch (here “patch” refers to the smallest unit of the repeating patterns). Minimum angle in a mesh is, in fact, one of the concerns of those who work on mesh generation part of an adaptive procedure.
-Figures 10 and 11 demonstrate the results obtained from application of SPR to heat and elasticity problems.  Interesting to note is that the results obtained from all patterns consisting of 2D Chevron pattern, at least at one face, lose superconvergent property.  Similar conclusion may be made from application of REP to the problem as shown in Figures 12 and 13.
+{| style="width: 100%;border-collapse: collapse;"
+|-
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_8233_Fig20a.png]]
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_9518_Fig20b.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(a)
+|  style="text-align: center;vertical-align: top;"|(b)
+|-
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_5602_Fig20c.png]]
+|  style="text-align: center;vertical-align: top;"|
+[[File:Draft_Samper_725107104_4543_Fig20d.png]]
+|-
+|  style="text-align: center;vertical-align: top;"|(c)
+|  style="text-align: center;vertical-align: top;"|(d)
+|}
-It may be noticed that existence of superconvergence is independent of the average of minimum angles since there are some cases in which the average minimum angle are quite small while the results are still superconvergent.   This means that, regarding error estimation part in an adaptive procedure, a mesh generation program should face more limitations with respect to mesh pattern.  On the other hand, if the results of recovery methods are taken as approximate gradients, the restrictions regarding minimum angle, i.e. skewness of elements, may be relaxed.
-Since we have so far studied structured mesh configurations, a question regarding the structure of the output mesh from a 3D mesh generation program, structured/unstructured mesh, may now be raised. Obviously, it is not possible to consider all possible unstructured mesh configurations and thus the answer needs an engineering judgment for generalization of the present results to unstructured meshes though the need of more studies, in this regard, is felt. Nevertheless, in this study, we examine the performance of the recovery methods in some artificially produced unstructured patterns.  For this we select the pattern of RRR and change the location of central node to produce skew elements inside the patch.   The results obtained, for some specific locations of the inside node with respect to the center of the patch are shown in Figure 14.  In this figure the patch is considered to be at interior parts of a fictitious domain. It can be seen that the robustness values are still in an acceptable range when the mesh becomes very skew though SPR exhibits better robustness for highly distorted meshes.
+Figure 20.  Convergence of the solution with Neumann and Dirichlet boundary condition-
+<math display="inline">u_1</math> and  <math display="inline">u_5</math> are nodal
+values of node one and five in Fig 5-c for the corner patch of the domain, and
+<math display="inline">\sigma </math> is the gradient at center of element
+-2-3-4; (a)/(c) Convergence of solution with respect to the number of Gauss points
+in a  <math display="inline">20\times 20</math> macro patch, (b)/(d) Convergence
+of solution with respect to the number of patches used along each direction while the
+number of Gauss points is fixed to 30 for all solutions.
-==5.  Conclusions==
+=5.  CONCLUSIONS=
-In this part of the paper we have presented a study on asymptotic behavior of two recovery methods, i.e. SPR and REP, in producing superconvergent gradients for three dimensional heat and elasticity problems.   It has been shown that application of both methods leads to similar robustness values either for interior of the domain or near boundaries.
+The formulation introduced by Babuška and co workers [13], for assessment of error estimator near flat boundaries, has been extended to kinked boundaries.   The new formulation given is suitable for study of either asymptotic behavior of a finite element solution near boundaries or general FEM solution for unbounded domains.
-Full superconvergent property is lost for tetrahedral elements when used near boundaries.  This effect is much pronounced for results from elasticity problems.  However, for brick elements the results exhibit the property in an asymptotic sense.  One of the reasons for loss of full superconvergent property seems to be the definitions of the recoveries near boundaries which basically use extrapolated information from polynomials fitted at inner patches.  Another reason is the fact that in both methods the information used for recovery procedure is extracted from the region in which boundary layer has not sufficiently decayed.   This effect seems to be much pronounced for tetrahedral elements since for brick elements asymptotic superconvergent property is still seen.  It may be concluded that some further research is needed to find better definitions for recovery methods near boundaries so that either the information from inner layers of recovered nodal values is utilized or the spurious boundary layer modes are captured and removed during the recovery procedures.  The method proposed in this paper provides a suitable tool for such studies.
+Asymptotic finite element solution consists of two main parts; one obtained from patches located at interior of the domain, and another obtained from a boundary layer solution irrespective of other boundaries at far ends of the domain.  The latter part of the solution decays towards the interior of the domain and, in this respect, is similar to a finite element solution of an unbounded domain.
-Some complementary results have been given for tests results for interior parts when various combinations of patters appear at three faces of a cuboid.  The conclusion is that all configurations produced by Chevron pattern at least in one face lose the full superconvergent effect.  This means that in mesh generation programs should face more restrictions regarding the error estimation part in adaptive solutions.  Some more results on artificially produced unstructured meshes have been given to examine the performance the recoveries in an asymptotic manner.  It is seen  that even for highly distorted mesh configurations the robustness indices remain in an acceptable range for both recovery methods considered.
+The proposed method is based on assuming proportionality relations between the finite element nodal values along two or three directions, when the mesh is of a repeatable pattern.  Writing the equilibrium equations over a unit patch leads to a characteristic equation containing all proportionality factors from which one may be expressed in terms of the others.  The boundary layer solution is obtained by integral of extended vectors arranged in spectral form from null space vectors of the equilibrium equations.  The decay condition is then applied by defining a feasible domain for proportionality factors.
-==References==
+Performance of the method has been examined for solution of unbounded domains in order to ensure its applicability to producing asymptotic finite element solution in an assessment procedure.  A sample problem has been solved using two forms of patterns to demonstrate the capability of the method in  reproducing the boundary layer solution. The convergence of the procedure has been studied with respect to number of patches as well as the number of integration points used.  The results show fast convergence for both cases.
-[1]  Boroomand, B. and Mossiaby, F., “On Application of robustness test for error estimators to patches near kinked boundaries in three dimensional problems: Part I, Formulation”, CIMNE, June 2004.
+In the next part of the paper we shall give the step by step procedure for robustness test followed by numerous results obtained from application of the test to some recovery based error estimators.
-[2]  Babuška, I. and Rheinboldt, W.C., “A posteriori error estimates for finite element method”, ''Int. J. Numer. Meth. Engrg.'', Vol. 12, pp. 1597-1615, 1978.
+==APPENDIX==
-[3]  Zienkiewicz, O.C. and Zhu, Z., “A simple error estimator and adaptivity procedure for practical engineering analysis”, ''Int. J. Numer. Meth. Engrg.'', Vol. 24, pp. 337-357, 1987.
+In this appendix we explain the reason for existence of proportionality property used in (52-a) to (52-c).
-[4]  Zienkiewicz, O.C. and Zhu, Z., “The superconvergent patch recovery and a posteriori error estimates. Part I: The recovery technique”, ''Int. J. Numer. Meth. Engrg.'', Vol. 33, pp. 1331-1364, 1992.
+Integral Equation (27) may be viewed as a finite element solution of a corresponding differential equation with constant coefficient as the following one
-[5]  Zienkiewicz, O.C. and Zhu, Z., “The superconvergent patch recovery and a posteriori error estimates. Part II: Error estimates and adaptivity”,  ''Int. J. Numer. Meth. Engrg.'', Vol. 33, pp. 1365-1382, 1992.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math>\boldsymbol{S}^T\boldsymbol{DS}\boldsymbol{u}_{h0}^p=</math><math>\boldsymbol{0}</math> on     <math>{\Omega }^{{'}}</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (A-1)
+|}
-[6]   Wiberg, N. E. and Abdulwahab, F., “Patch recovery based on superconvergent derivatives and equilibrium”, ''Int. J. Numer. Meth. Engrg.'', Vol. 36, pp. 2703-2724, 1993.
-[7]    Blacker, T.D. and Belytschko, T., “Superconvergent patch recovery with equilibrium and conjoint interpolation enhancement’ ''Int. J. Numer. Meth. Engrg.'', Vol. 37, pp. 517-536, 1994.
+For scalar problems such as heat conduction ones, for instance, the differential equation becomes
-[8] Wiberg, N. E., Abdulwahab, F. and Ziukas, S., “Enhanced superconvergent patch recovery incorporation equilibrium and boundary condition”, ''Int. J. Numer. Meth. Engrg.'', Vol. 37, pp. 3417-3440, 1994.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math>{\nabla }^2u=0</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (A-2)
+|}
-[9]  Babuška, I., Strouboulis, T. and Upadhyay, C.S., “ A model study of the quality of a posteriori error estimaors for linear elliptic problems. Error estimation interior of patchwise uniform grids of triangles”, ''Comput. Methods Appl. Mech. Eng.'', Vol. 114, pp. 307-378, 1994.
+forgetting the subscript and superscripts for simplicity.  Since the equation is of constant coefficient type, the solution may be expressed in an exponential form as  <math display="inline">u=ce^{\alpha x+\beta y}</math> .  This leads to a characteristic equation as   <math display="inline">{\alpha }^2+{\beta }^2=0</math> after substitution in the equation.  The exponential form of solution of (A-2) may be written as (see also [17])
-[10]  Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K., “Validation of a posteriori error estimators by numerical approach”, ''Int. J. Numer. Meth. Engrg.'', Vol. 37, pp. 1073-1123, 1994.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+| <math>u={\int }_{\mbox{ }\alpha }c_{\left(\alpha \right)}e^{\alpha x+\beta y}d\alpha \qquad \beta =f(\alpha )</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (A-3)
+|}
-[11]  Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K. “An objective criterion for assessing the reliability of a posteriori error estimators in finite element computations”, ''U.S.A.C.M. Bulletin'', No.7 , pp. 4-16, 1994.
+or
-[12]   Babuška, I., Strouboulis, T., Mathur, A. and Upadhyay, C.S., “Pollution-error in the h-version of the finite element method and the local quality of a posteriori error estimators”, ''Finite Elem. Anal. Des., ''Vol. 17, pp. 237-321, 1994.
+{| class='formulaSCP' style='width: 100%;'
+|-
+|
+{| style='margin:auto;width: 100%; text-align:center;'
+|-
+|  <math display="inline">u={\int }_{\mbox{ }\beta }c_{\left(\beta \right)}e^{\alpha x+\beta y}d\beta </math>
+| <math display="inline">\alpha =g(\beta )</math>
+|}
+| style='width: 5px;text-align: right;white-space: nowrap;' | (A-4)
+|}
-[13]  Babuška, I., Strouboulis, T., Upadhyay, C.S. and Gangaraj, S.K., “A posteriori estimation and adaptive control of the pollution error in the h-version of the finite element method”,  ''Int. J. Numer. Meth. Engrg.'', Vol. 38, pp. 4207-4235, 1995.
+where integrations are taken on feasible domain of  <math display="inline">\alpha </math> or  <math display="inline">\beta </math> , e.g. conditions for decay property. Considering single solution from (A-3) for example, it may be seen that the following relations exist
-[14]  Babuška, I., Strouboulis, T., Gangaraj, S.K. and Upadhyay, C.S., “Pollution error in the h-version of the finite element method and the local quality of the recovered derivatives”,  ''Comput. Methods Appl. Mech. Eng.'', Vol. 140, pp. 1-39, 1997.
-[15]  Boroomand, B. and Zienkiewicz, O.C., “Recovery by equilibrium in patches (REP)”, ''Int. J. Numer. Meth. Engrg.'', Vol. 40, pp. 137-164, 1997.
+{|
+|-
+| <math>u_{\left(x,y\right)}=c_{\left(\alpha \right)}e^{\alpha x+\beta y}</math>
+| <math>u_{\left(x,y+L_y\right)}=c_{\left(\alpha \right)}e^{\alpha x+\beta (y+L_y)}=</math><math>c_{\left(\alpha \right)}e^{(\alpha x+\beta y)+\beta L_y}=</math><math>c_{\left(\alpha \right)}e^{\alpha x+\beta y}e^{\beta L_y}=</math><math>{\mu }_{02}u_{\left(x,y\right)}</math>
+|}
-[16]  Boroomand, B. and Zienkiewicz, O.C., “An improved REP recovery and efectivity robustness test”, ''Int. J. Numer. Meth. Engrg.'', Vol. 33, pp. 3247-3277, 1997.
+{|
+|-
+| <math>u_{\left(x,y+mL_y\right)}=c_{\left(\alpha \right)}e^{\alpha x+\beta (y+mL_y)}=</math><math>c_{\left(\alpha \right)}e^{(\alpha x+\beta y)+m\beta L_y}=</math><math>c_{\left(\alpha \right)}e^{\alpha x+\beta y}e^{m\beta L_y}=</math><math>{\left({\mu }_{02}\right)}^mu_{\left(x,y\right)}</math>
+|}
-[17]  Zienkiewicz,  O.C., Boroomand, B. and Zhu, Z., “Recovery procedures in error estimation and adaptivity; Part I: Adaptivity in linear problems”,'' Comput. Methods Appl. Mech. Eng.'', Vol. 176, pp. 111-125, 1999.
+Similarly
-[18]  Boroomand, B. and Zienkiewicz, O.C., “Recovery procedures in error estimation and adaptivity; Part II: Adaptivity in nonlinear problrms of elaso-plasticity behavior”, ''Comput. Methods Appl. Mech. Eng.'', Vol. 176, pp. 127-16, 1999.
-[19]  Lo, S.H. and Lee, C.K., “On using different recovery procedures for the construction of smoothed stress in finite element method”, ''Int. J. Numer. Meth. Engrg.'', Vol. 43, pp. 1223-1252, 1998.
+{|
+|-
+| <math>u_{\left(x+L_x,y\right)}=c_{\left(\alpha \right)}e^{\alpha (x+L_x)+\beta y}=</math><math>{\mu }_{01}u_{\left(x,y\right)}</math>
+| <math>u_{\left(x+nL_x,y\right)}=c_{\left(\alpha \right)}e^{\alpha (x+nL_x)+\beta y}=</math><math>{\left({\mu }_{01}\right)}^nu_{\left(x,y\right)}</math>
+|}
-[20]  Babuška, I., Strouboulis, T. and Upadhyay, C.S., “A model study of the quality of a posteriori error estimators for finite element solutions of linear elliptic problems, with particular reference to the behavior near the boundary”,  ''Int. J. Numer. Meth. Engrg.'', Vol. 40, pp. 2521-2577, 1997.
+or more generally
-[21]  Zienkiewicz, O.C. and Zhu, Z., “The superconvergent patch recovery (SPR) and adaptive finite element refinement”, ''Comput. Methods Appl. Mech. Eng.'', Vol. 101, pp. 207-224, 1992.
+<math>u_{\left(x+nL_x,y+mL_y\right)}=c_{\left(\alpha \right)}e^{\alpha (x+nL_x)+\beta (y+mL_y)}=</math><math>{\left({\mu }_{01}\right)}^n{\left({\mu }_{02}\right)}^mu_{\left(x,y\right)}</math>
-[22]  Zienkiewicz, O.C., Zhu, Z. and Wu, J., “Superconvergent recovery technique-some further tests”’ ''Commun. Num. Meth. Eng.'', Vol. 9, pp.251-258,1993.
-[23]  Zienkiewicz, O.C. and Zhu, Z., “Superconvergence and the superconvergent patch recovery”, ''Finite Elem. Anal. Design.'', Vol. 19, pp.11-23, 1995.
+Which is consistent with assumption made in (52-a) to (52-c).  We note that here the differential equation is satisfied exactly whereas in relations (52) the aim is to satisfy the integral form of the equation.   Of course the proportionality factors in the two approaches are not similar since one is an approximate finite element solution of another.
+An important observation may be made by looking in to the characteristic equation  <math display="inline">{\alpha }^2+{\beta }^2=0</math> . It can be seen that for a given value of  <math display="inline">\alpha </math> there will be two values for  <math display="inline">\beta </math> one of which corresponds to decay condition, i.e.  <math display="inline">e^{-i\alpha y}</math> or  <math display="inline">e^{i\alpha y}</math> (note that  <math display="inline">\alpha </math> can be a complex value).  Hence for a given value of  <math display="inline">{\mu }_{01}</math> , two values as  <math display="inline">{\mu }_{02}</math> and  <math display="inline">{\mu }_{02}^{-1}</math> are both feasible but the one with absolute value less than unity corresponds to decay condition. Similar relations may be written for a general three dimensional problem.
+==References==
+[1]  Oden, J.T., Belytschko, T., Babuška, I. and Hughes, T.J.R., “Research directions in computational mechanics” ''Comput. Methods Appl. Mech. Eng.'', Vol. 192, pp. 913-922, 2003.
+[2]  Zienkiewicz, O.C., “Achievements and some unsolved problems of finite element method” , ''Int. J. Numer. Meth. Engrg.'', Vol. 47, pp. 9-28, 2000.
+[3]  Hinton, E. and Campbell, J.S., “Local and global smoothing of discontinuous finite element functions using a least square method”, ''Int. J. Numer. Meth. Engrg.'', Vol. 8, pp. 461-480, 1976.
+[4]  Zienkiewicz, O.C. and Zhu, Z., “A simple error estimator and adaptivity procedure for practical engineering analysis”, ''Int. J. Numer. Meth. Engrg.'', Vol. 24, pp. 337-357, 1987.
+[5]  Tessler, A., Riggs, H.R. and Dambach, M., “A novel four-node quadrilateral smoothing element for stress enhancement and error estimation”, ''Int. J. Numer. Meth. Engrg.'', Vol. 44, pp. 1527-1543, 1999.
+[6]  Hiller, J.F. and Bathe, K.J., “On higher-order-accuracy points in isotropic finite element analysis and application to error assessment”, ''Comput. & Struc.,'' Vol. 79, pp. 1275-1285, 2001.
+[7]  Zienkiewicz, O.C. and Zhu, Z., “The superconvergent patch recovery and a posteriori error estimates. Part II: Error estimates and adaptivity”,  ''Int. J. Numer. Meth. Engrg.'', Vol. 33, pp. 1365-1382, 1992.
+[8]  Babuška, I., Strouboulis, T. and Upadhyay, C.S., “ A model study of the quality of a posteriori error estimaors for linear elliptic problems. Error estimation interior of patchwise uniform grids of triangles”, ''Comput. Methods Appl. Mech. Eng.'', Vol. 114, pp. 307-378, 1994.
+[9]  Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K., “Validation of a posteriori error estimators by numerical approach”, ''Int. J. Numer. Meth. Engrg.'', Vol. 37, pp. 1073-1123, 1994.
+[10]    Babuška, I., Strouboulis, T., Upadhyay, C.S., Gangaraj, S.K. and Copps, K. “An objective criterion for assessing the reliability of a posteriori error estimators in finite element computations”, ''U.S.A.C.M. Bulletin'', No.7 , pp. 4-16, 1994.
+[11] Zhang, Z. and Zhang, S., “ Derivative superconvergence of rectangular finite elements for the Reissner-Mindlin plate”, ''Comput. Methods Appl. Mech. Eng., ''Vol.  134, pp. 1-16, 1996.
+[12] Zhang, Z., “ Superconvergence in projected –shear plate-bending finite element methods”, ''Numer. Meth. Partial Diff. Eqn., ''Vol.  14, pp. 367-386, 1998.
+[13]  Babuška, I., Strouboulis, T. and Upadhyay, C.S., “A model study of the quality of a posteriori error estimators for finite element solutions of linear elliptic problems, with particular reference to the behavior near the boundary”,  ''Int. J. Numer. Meth. Engrg.'', Vol. 40, pp. 2521-2577, 1997.
+[14]  Ainsworth, M. and Oden, J.T., ''A posteriori error estimation in finite element analysis'', John Wiley, 2000.
+[15]  Zienkiewicz, O. C. and Taylor, R. L., “''The finite element method''”, 5<sup>th</sup> edition, Butterworth-Heineman, 2000.
+[16]  Boroomand, B. and Zienkiewicz, O.C., “An improved REP recovery and efectivity robustness test”, ''Int. J. Numer. Meth. Engrg.'', Vol. 33, pp. 3247-3277, 1997.
-[24]  Zienkiewicz, O. C. and Taylor, R. L., “''The finite element method''”, 5<sup>th</sup> edition, Butterworth-Heineman, 2000.
+[17] Hildebrand, F., ''Advanced calculus for applications,'' Second edition, Prentice-Hall, 1976.