Abstract:
The present invention relates to an estimation of computational errors appearing in the finite element calculations, particularly to method, system and program product for verification of the accuracy of numerical data measured in terms of problem-oriented criteria, where the numerical data are computed in the process of solution of a boundary value problem. In addition to the primal problem a certain adjoin problem is formed and solved. The method is based on two principles:  
     (a) the original and adjoint problems are solved on non-coinciding meshes, and  
     (b) the term presenting the product of errors arising in the primal and adjoint problems is estimated by the gradient recovery technique.

Description:
BACKGROUND OF THE INVENTION  
         [0001]    1. Field of the Invention  
           [0002]    The present invention relates to an estimation of computational errors appearing in the finite element calculations, particularly to method, system and program product for verification of the accuracy of numerical data measured in terms of problem-oriented criteria, where the numerical data are computed in the process of solution of a boundary value problem, which is defined by one or several partial differential equations (PDEs), governing an unknown physical quantity in a solution domain, and boundary conditions.  
           [0003]    2. Description of the Prior Art  
           [0004]    Finite Element Method (FEM) is a very powerful and nowadays the most popular numerical method for solving the partial differential equations (PDEs). The concept of FEM consists of discretizing the solution domain into a set (called the mesh) consisting of small elements with simple shapes. The mesh is characterized by a number of degrees of freedom (DOF), we denote it by n in the text. In structural mechanics, the DOF represent the virtual displacements on points (nodes), where the (finite) elements are constructed. The finite element (FE) solution is built from the solution of a system of n linear algebraic equations. Engineers have years of experience and practice in the finite element calculations and analysis. When designing e.g. a new structure, the mechanical engineer can use the FEM to calculate a tension distribution in a structure (e.g., in frame, crankshaft or beam). Normally, there are a few critical areas to be carefully analysed. An ordinary approach is to perform computations on more and more dense mesh, hoping that the computed approximate solutions are becoming more and more closer to the exact solution (which is normally not known at all). But the closeness of such calculated approximations to the exact solution can be illusory in many respects, and the computed data (numbers) can be far from the truth, thus being rather unreliable and useless. A solution of this problem is an availability of certain tools (software) for the engineers and designers for explicit control of the error between the exact solution and the computed approximation. A development of theoretical backgrounds and various techniques for building of such tools (called a priori error estimates) has been a main topic in the modern numerical analysis last two decades. For the finite element methods, such a priori error estimates are usually obtained either by estimating a weak norm of the residual [1] or by using special postproces sing procedures [5]. For the Galerkin (finite element) approximations of linear elliptic problems, they estimate the error in the global (energy) norm and also provide an error indicators that are further used in various mesh adaptive procedures. Global error estimates give a general presentation on the quality of an approximate solution and a stopping criteria. However, from the viewpoint of engineering purposes, such an information is often not sufficient. In many cases, the engineers are interested not in the value of the total error, but, as was mentioned above, in errors over certain critical areas: subdomains, lines, or at special points. A possible way of estimating such type of errors is to introduce a certain linear functional l associated with a critical area, where l(u) is often called as a problem-oriented criterion or “quantity of interest”, and to obtain an estimate for the value of l(u−v), where u is the exact solution and v is the approximate one. Known methods [2] find estimates for l(u−u h ) for the Galerkin approximations u h  by employing an auxiliary (often called as adjoint) problem, whose right-hand side is formed by the functional l. If the Galerkin approximation of the adjoint problem is computed on the same mesh as u h , then the functional value l(u−u h ) is expressed via a certain integral functional, which can be estimated by using, e.g., the “equilibrated residual method” [4]. However, all the results obtained by now are rather far from being satisfactory.  
           [0005]    What is needed is a reliable and effective method, system and program product for verification of the accuracy of numerical data measured in terms of problem-oriented criteria.  
         OBJECTS AND SUMMARY OF THE INVENTION  
         [0006]    We present a new method for estimation of the “quantities of interest” [3] as well as a system and a program product for implementing the method. It is based on two principles:  
           [0007]    (a) the original and adjoint problems are solved on non-coinciding meshes, and  
           [0008]    (b) the term presenting the product of errors arising in the primal and adjoint problems is estimated by the gradient recovery technique widely used in various applied problems.  
           [0009]    These features differ our method from others, who are usually assumed that the Galerkin approximations of the primal and auxiliary (adjoint) problems should be computed on the same finite dimensional subspaces.  
           [0010]    Using the fact of non-coincidence of the meshes, we obtain a new estimator which is a sum of two terms (E 0 ) and (E 1 ), where the first term (E 0 ) contains the major part of the “quantity of interest”, but vanishes whenever the mesh coincide. The effectivity of the method is demonstrated for one- and two-dimensional problems, part of which is presented in the text later on.  
           [0011]    In the preferred embodiment the estimator (E) is a sum of the first term (E 0 ) and the second term (E 1 ), where  
           [0012]    the first term (E 0 ) is computable as a sum of three parts, where the first part is an integral over the solution domain (Ω) of a product of the source function (f) and the second approximation (v τ ), the second part is an integral over the Neumann boundary (Γ 2 ) of a product of the boundary function (g) and the second approximation (v τ ), and the third part is a minus integral over the solution domain (Ω) of an energy scalar products of the gradients of the first (u h ) and the second (v τ ) approximations, and  
           [0013]    the second term (E 1 ) is the integral over the solution domain (Ω) of the energy scalar product of two vector-functions, where the first vector-function is the difference between the gradient of the first approximation (u h ) and the averaged gradient of the first approximation (u h ), and the second vector-function is the difference between the gradient of the second approximation (v τ ) and the averaged gradient of the second approximation (v τ ).  
           [0014]    The effectivity of the method proposed, strongly increases when one is interested not in a single solution of the primal problem for a concrete data, but analyzes a series of approximate solutions for a certain set of boundary conditions and various right-hand sides (and such a situation is very typical in the engineering design, when it is necessary to model the behaviour of a construction for various working regimes). In this case, the adjoint problem must be solved only once for each “quantity of interest”, and its solution can be further used in testing the accuracy of approximate solutions of various primal problems (i.e. with various source functions).  
           [0015]    The present invention can be implemented using a computer. An exemplary computer includes one or more processors, memory (like ROM, RAM and hard-drive), input device (like keyboard, floppy disk drive) and output device (like display and printer). Method can be implemented directly for any software using FEM as a computational tool. Without any new software the first and second approximations as well as computing the estimator must be executed separately. We present new software, which includes all necessary routines and interface for an implementation. 
       
    
    
     BRIEF DESCRIPTION OF THE DRAWINGS  
       [0016]    In the following, the invention is examined with the aid of the accompanying drawings and examples. In the drawings illustrates a solution domain with basic denotations for the diffusion type problem  
         [0017]    [0017]FIG. 1 illustrates a solution domain with basic denotations  
         [0018]    [0018]FIG. 2 illustrates a part of standard finite element mesh and a patch associated with one of its nodes  
         [0019]    [0019]FIG. 3 illustrates an elastic body as a solution domain with basic denotations for the linear elasticity problem  
         [0020]    [0020]FIG. 4 illustrates main desktop of CHECKER-software  
         [0021]    [0021]FIG. 5 illustrates setting the domain and the zone of interest in Matlab® PDE Toolbox  
         [0022]    [0022]FIG. 6 illustrates the export window of Matlab® PDE Toolbox for exporting domain data and formula data for the primal problem to CHECKER-software  
         [0023]    [0023]FIG. 7 illustrates setting the boundary conditions in Matlab® PDE Toolbox  
         [0024]    [0024]FIG. 8 illustrates the export window for exporting data of the zone of interest and boundary conditions to CHECKER-software  
         [0025]    [0025]FIG. 9 illustrates the process of solving the primal problem in CHECKER-software  
         [0026]    [0026]FIG. 10 illustrates the process of solving the auxiliary (adjoint) problem in CHECKER-software  
         [0027]    [0027]FIG. 11 illustrates the process of solving the auxiliary problem with a local refinement in the zone of interest in CHECKER-software  
         [0028]    [0028]FIG. 12 illustrates CHECKER desktop with visible results of  
         [0029]    [0029]FIG. 13 illustrates a block diagram of CHECKER-software with short step references A-E  
         [0030]    [0030]FIG. 14 illustrates the source function of the one-dimensional test problem  
         [0031]    [0031]FIG. 15 illustrates the exact and the approximate solutions and the zone of interest of the one-dimensional test problem  
         [0032]    [0032]FIG. 16 illustrates the behaviour of the effectivity index in various mesh parameters  
         [0033]    [0033]FIG. 17 illustrates the meshes T h  and T τ  for the one-dimensional test problem  
         [0034]    [0034]FIG. 18 illustrates a flowchart for computing an approximation and an error estimation 
     
    
     DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS  
       [0035]    Table of Contents  
         [0036]    I. Model Theory  
         [0037]    A. Introduction  
         [0038]    B. Problem Formulation—Diffusion-Type Problem  
         [0039]    1. Mathematical model  
         [0040]    2. Finite element solution  
         [0041]    3. Problem-oriented criterion  
         [0042]    4. Error estimation  
         [0043]    5. Technique for the error estimation in terms of problem-oriented criteria  
         [0044]    6. Auxiliary problem and its finite element solution  
         [0045]    7. Gradient averaging procedures  
         [0046]    8. Formula for calculation of the estimator  
         [0047]    9. Error estimation algorithm  
         [0048]    C. Problem Formulation—Linear Elasticity Problem  
         [0049]    1. Mathematical model  
         [0050]    2. Finite element solution  
         [0051]    3. Problem-oriented criterion  
         [0052]    4. Error estimation  
         [0053]    5. Technique for the error estimation in terms of problem-oriented criteria  
         [0054]    6. Auxiliary problem and its finite element solution  
         [0055]    7. Tensor averaging procedures  
         [0056]    8. Formula for calculation of the estimator  
         [0057]    9. Error estimation algorithm  
         [0058]    D. Advantages Of The Technique  
         [0059]    II. Implementation In A System  
         [0060]    A. Introduction  
         [0061]    B. Presentation of Software Product with a Model  
         [0062]    III. Model Example—One Dimensional Case  
         [0063]    I Model Theory  
       A—Introduction  
       [0064]    The method presented below is designed to verify the accuracy of FE approximations, measured in terms of “problem-oriented” criterion that can be chosen by users. The method is applic able to various problems in mechanics and physics embracing diffusion and linear elasticity problems. It can be easily coded and attached as an independent programme-checker to the most of existing educational and industrial codes.  
       B—Problem Formulation Part—Diffusion-Type Problem  
       [0065]    B.1 Mathematical Model  
         [0066]    We consider physical process that can be formally modelled in the form of a boundary value problem of elliptic type (BVP) as follows referring to FIG. 1: Find a function u=u(x 1 , . . . ,x d ) of d variables x 1 , . . . ,x d  such that  
                 -       ∑     i   ,     j   =   1       d              ∂               ∂     x   i              (       a   ij            ∂   u       ∂     x   j           )           =     f                 in                 Ω       ,           (   1   )                               
 u=u 0  on Γ 1 , (2a)  
                 ∑     i   ,     j   =   1       d            a   ij            ∂   u       ∂     x   j              n   i         =     g                 on                     Γ   2     .               (     2      b     )                               
 
         [0067]    In the above, relation (1) is a process governing equation in the solution domain Ω⊂R d , where a ij =a ij (x 1 , . . . ,x d ),i,j=1, . . . ,d, are the given coefficients of the problem that usually describe the diffusion properties of the respective media. This model also applies to stationary magnetic, electric, and temperature field, or some phenomena in the linear elasticity or fluid flow. The function f=f(x 1 , . . . ,x d ) can be viewed as a given source function. The coefficients are assumed to be bounded functions and function f is assumed to be square summable. Further, relations (2a) and (2b) are boundary conditions on the boundary Γ of Ω, more precisely, they prescribe a behaviour of the solution and its derivatives on two nonintersecting boundary parts Γ 1  and Γ 2 , respectively, i.e., Γ=Γ 1∪Γ   2 . The function u 0 =u 0 (x 1 , . . . ,x d ) belongs to the natural energy class of the problem and g=g(x 1 , . . . ,x d ) is a square summable function, symbol n i  denotes the i-th component of the outward normal vector n=n(x 1 , . . . ,x d ) to the boundary Γ, i.e., n=[n 1 , . . . ,n d ] T  (see FIG. 1).  
         [0068]    The above assumptions are not restrictive and cover the bulk of practically meaningful cases.  
         [0069]    B.2 Finite Element Solution of (1)-(2)  
         [0070]    Let V h  be a finite-dimensional space constructed by means of a selected set of finite element trial functions defined on a standard finite element mesh T h  (called the first mesh) over the solution domain Ω, see FIG. 2. We notice that space V h  is chosen so that its functions w h  vanish on Γ 1 .  
         [0071]    The first approximation for problem (1)-(2) is defined as a function  
           u   h   =u   h ( x   1   , . . . ,x   d )ε V   h   +u   0 ,  
         [0072]    such that  
                 ∫   Ω            ∑     i   ,     j   =   1       d            a   ij            ∂     u   h         ∂     x   j                ∂     w   h         ∂     x   i                 x           =         ∫   Ω            fw   h             x         +       ∫     Γ   2              gw   h             s                     ∀       w   h     ∈       V   h     .                       (   3   )                               
 
         [0073]    B.3 Problem-Oriented Criterion  
         [0074]    Users are often interested not only in the overall error e=u−u h , but also in its local behaviour, e.g., in a certain subdomain ω⊂Ω. One way to get information about the local behaviour of u−u h  is to measure the error in terms of specially selected problem-oriented criteria on “quantity of interest”. One of the most typical criteria of such a type is presented by the integral  
                 ∫   Ω            ϕ        (     u   -     u   h       )               x         ,           (   4   )                               
 
         [0075]    where φ=φ(x 1 . . . ,x d ) is a selected weight function vanishing outside of the zone of interest.  
         [0076]    B.4 Error Estimation  
         [0077]    Error estimators are intended to show the value of the error in the subdomain, which is the most interesting for an engineer. Also they suggest information on that how to improve an approximate solution computed if the accuracy obtained is not sufficient. At this point we note that our estimator is given by the integral whose integrand could serve as the respective indicator. Namely, if the accuracy achieved is not sufficient, then additional degrees of freedom (new elements) must be added in those parts of the solution domain, where the integrand is excessively high.  
         [0078]    B.5 Technique for the Error Estimation in Terms of Problem Oriented Criteria  
         [0079]    To present the technique for estimation of criterion (4), we describe two technical problems that must be previously solved. They consist of finding an approximate solution of an auxiliary problem (so-called adjoint problem) and making a certain post-processing of this solution v o , and also making a post-processing of the finite element approximation u h .  
         [0080]    B.6 Auxiliary Problem and its Finite Element Solution  
         [0081]    Let V τ  be the second finite-dimensional space constructed by means of a selected set of finite element trial functions on another standard finite element mesh T τ  (called the second mesh) over the solution domain Ω. We notice that space V τ  is chosen so that its functions w r  vanish on Γ, and also that T τ  need not to coincide with T h , see FIG. 2. A part of standard finite element mesh and a patch associated with node.  
         [0082]    Consider the auxiliary finite-dimensional problem as follows: Find a function  
           v   τ   =v   τ ( x   1   , . . . , x   d )δ V   τ ,  
         [0083]    such that  
                 ∫   Ω            ∑     i   ,     j   =   1       d            a   ji            ∂     v   τ         ∂     x   j                ∂     w   τ         ∂     x   i                 x           =       ∫   Ω          ϕ                   w   τ             x                     ∀       w   τ     ∈       V   τ     .                     (   5   )                               
 
         [0084]    B.7 Gradient Averaging Procedures  
         [0085]    On the first mesh T h , we define the gradient averaging transformation G h  mapping the gradient of the finite element approximation u h  
                 ∇     u   h       =       [         ∂     u   h         ∂     x   1         ,   …              ,       ∂     u   h         ∂     x   d           ]     T       ,           (   6   )                               
 
         [0086]    which is constant over each element of the finite element mesh, into a vector-valued continuous piecewise affine function  
           G   h (∇u h )=[ G   h   1 (∇ u   h ), . . . ,  G   h   d (∇u h )] T ,  (7)  
         [0087]    by setting each its nodal value as the mean value of ∇u h  on all elements of the patch P(x o ) associated with corresponding node x o  in the mesh T h .  
         [0088]    More precisely, let  
           ∂     u   h         ∂     x   i                   T   k                             
 
         [0089]    denote the value of i-th coordinate of gradient ∇u h  over triangle T k , let x o  be one of nodes of finite element mesh T h , and let T 1 , . . . , T N     xo    be elements of the mesh having node x o  as one of their vertices. The union of such elements is called the patch and denoted as P(x o ). (Thus, in terms of FIG. 2, the patch consists of six triangles, i.e., N x =6.) Then, we define  
                   G   h   i          (     ∇     u   h       )            (     x   *     )       =       1     N     x   *                ∑       T   k     ∈     P        (     x   *     )                            ∂     u   h         ∂     x   i                     T   k            ,     i   =   1     ,   …              ,   d   ,     
        or                     (     8      a     )                     G   h   i          (     ∇     u   h       )            (     x   *     )       =       1       measT   1     +   …   +     measT     N     x   *                    ∑       T   k     ∈     P        (     x   *     )                            measT   k     ·       ∂     u   h         ∂     x   i                       T   k            ,     i   =   1     ,   …              ,   d   ,                     (     8      b     )                               
 
         [0090]    where the symbol meas denotes the area of corresponding triangle.  
         [0091]    Having G h (∇u h ) defined at all nodes of the mesh T h , we uniquely define continuous piecewise affine function G h   i (∇u h ) over the whole domain Ω. In this way, the vector-valued continuous piecewise affine function G h (∇u h ) in (7) is built. Similarly, on T τ , we define the gradient averaging transformation G τ  mapping the gradient of the finite element approximation V τ  
                 ∇     v   τ       =       [         ∂     v   τ         ∂     x   1         ,   …              ,       ∂     v   τ         ∂     x   d           ]     T       ,           (   9   )                               
 
         [0092]    which is constant over each element of the finite element mesh, into a vector-valued continuous piecewise affine function  
           G   τ (∇ v   τ )=[ G   τ   1 (∇ v   τ ), . . . , G   τ   d (∇ v   τ )] T ,  (10)  
         [0093]    by setting each its nodal value as the mean value of ∇v τ  on all elements of the patch P(x o ) associated with corresponding node x o  in the mesh T τ .  
         [0094]    More precisely, let  
         [0095]    denote the value of i-th coordinate of gradient ∇v τ  over triangle T k , let x o  be one of nodes of finite element mesh T τ , and let T 1 , . . . ,  
                   E   0          (       u   h     ,     v   τ       )       =         ∫   Ω          f                   v   τ             x         +       ∫     Γ   2              gv   τ             s         -       ∫   Ω            ∑     i   ,     j   =   1       d            a   ij            ∂     u   h         ∂     x   j                ∂     v   τ         ∂     x   i                 x               ,     
        and           (   13   )                   E   1          (       u   h     ,     v   τ       )       =       ∫   Ω            ∑     i   ,     j   =   1       d              a   ij          (         ∂     u   h         ∂     x   j         -       G   h   j          (     ∇     u   h       )         )            (         ∂     v   τ         ∂     x   i         -       G   τ   i          (     ∇     v   τ       )         )               x     .                   (   14   )                               
 
         [0096]    T N     xo    be elements of the mesh having node x o  as one of their vertices. The respective patch is denoted as P(x o ). Then  
                   G   τ   i          (     ∇     v   τ       )            (     x   o     )       =       1     N     x   o                ∑       T   k     ∈     P        (     x   o     )                            ∂     v   τ         ∂     x   i                     T   k            ,     i   =   1     ,   …              ,     d   .     
        or                       (     11      a     )                     G   τ   i          (     ∇     v   τ       )            (     x   o     )       =       1       measT   1     +   …   +     measT     N     x   o                    ∑       T   k     ∈     P        (     x   o     )                            measT   k     ·       ∂     v   τ         ∂     x   i                       T   k            ,     i   =   1     ,   …              ,     d   .                       (     11      b     )                               
 
         [0097]    Having G τ   i (∇v τ )(x o ) at all nodes x o , we uniquely define continuous piecewise affine function G τ   i (∇v τ ) over the whole domain Ω. In this way, the vector-valued continuous piecewise affine function G τ (∇v τ ) in (10) is built.  
         [0098]    Such averaging transformations are widely used in the finite element calculations. Usually they lead to computationally inexpensive algorithms.  
         [0099]    B.8 Formula for Calculation of the Estimator  
         [0100]    Our method estimates the quantity (4) by the quantity E(u h ,v τ ) presented by the following formula:  
           E ( u   h   ,v   τ ):= E   0 ( u   h   ,v   τ )+ E   1 ( u   h ,v τ ),  (12)  
         [0101]    where  
           ∂     v   τ         ∂     x   i                   T   k                             
 
         [0102]    The functional E(u h ,v τ ) is directly computable once the approximations uh and v, are found.  
         [0103]    B.9. Error Estimation Algorithm  
         [0104]    Step 1: Pose a problem of type (1)-(2a,b); (Step A in FIGS. 13 and 18)  
         [0105]    Step 2: Find u h  in (3); (Step B)  
         [0106]    Step 3: Select criterion (4), i.e. subdomain ω and function φ; (Step C)  
         [0107]    Step 4: Solve auxiliary problem (5) and find v τ ; (Step D)  
         [0108]    Step 5: Construct/form G h (∇u h ) by (8a) or (8b); (Estimation  
         [0109]    Step 6: Construct/form G τ (∇v τ ) by (11a) or (11b); (Estimation)  
         [0110]    Step 7: Compute E 0 (u h ,v τ ); (Estimation)  
         [0111]    Step 8: Compute E 1 (u h ,v τ ); (Estimation)  
         [0112]    Step 9: Compute E(u h ,v τ ); (Estimation)  
         [0113]    (Steps 5-9 correspond step E in FIGS. 13 and 18.  
         [0114]    Optional step 10: change initial parameters. (Step F1 or F2)  
       C—Problem Formulation—Linear Elasticity Problem  
       [0115]    C.1 Mathematical Model  
         [0116]    We consider mechanical problem of elasticity that can be formally presented in the form of a system of linear elasticity equations as follows: Find a vector-function (called the displacement) u=u(x 1 , . . . ,x d )=[u 1 , . . . ,u d ] T  of d variables x 1 , . . . ,x d  such that, for i=1, . . . ,d  
                 -       ∑     j   ,   k   ,     l   =   1       d              ∂               ∂     x   j              (       L   ijkl            ɛ   kl          (   u   )         )           =       f   i                   in                 Ω       ,           (   1   )                               
 u i =u 0   i  on Γ 1 ,  (2a)  
                 ∑     j   ,   k   ,     l   =   1       d            n   j          L   ijkl            ɛ   kl          (   u   )           =       g   i                   on                     Γ   2     .               (     2      b     )                               
 
         [0117]    In the above, relation (1) is a governing system of equations in the solution domain (body) Ω⊂R d , where L ijkl =L ijkl (x 1 , . . . ,x d ),i,j=1, . . . ,d, are the given coefficients of the problem that usually describe the elastic properties of the respective body, and the vector-function f=f(x 1 , . . . ,x d )=[f 1 , . . . f d ] T  represents a given body force. The coefficients are assumed to be bounded functions and function f is assumed to be square summable. In addition, the coefficients satisfy the following symmetry conditions  
         L jikl =L ijkl =L klij , i,j,k,l=1, . . . ,d.  (3)  
         [0118]    Further, the functions ε ij (u),i,j=1, . . . ,d, are defined as follows  
                   ɛ   ij          (   u   )       =       1   2          (         ∂     u   i         ∂     x   j         +       ∂     u   j         ∂     x   i           )         ,           (   4   )                               
 
         [0119]    they form the so-called small strain (or deformation) tensor (d×d matrix) ε(u).  
         [0120]    The relations (2a) and (2b) are boundary conditions on the boundary Γ of Ω, more precisely, they prescribe a behaviour of the solution and its derivatives on two nonintersecting boundary parts Γ 1  and Γ 2 , respectively, i.e., Γ=Γ 1 ∪Γ 2 . The vector-functions u 0 =u 0 (x 1 , . . . ,x d )=[u 0   1 , . . . ,u 0   d ] T  represents the external force and belongs to natural energy class of the problem and g=g(x 1 , . . . ,x d )=[g 1 , . . . ,g d ] T  is a square summable function. The function g represents the external force. The symbol n i  denotes the i-th component of the outward normal vector n=n(x 1 , . . . ,x d ) to the boundary Γ, i.e., n=[n 1 , . . . ,n d ] T  (see FIG. 3).  
         [0121]    Above assumptions are not restrictive and cover the bulk of practically meaningful cases.  
         [0122]    C.2 Finite Element Solution of (1)-(2)  
         [0123]    Let V h  be a finite-dimensional space constructed by means of a selected set of finite element trial (d-dimensional) vector-functions defined on commonly-used finite element mesh T h  over Ω. We notice that space V h  is chosen so that its functions w h  vanish on Γ 1 .  
         [0124]    The finite element approximation for problem (1)-(2) is defined then as a vector-function such that  
           u   h   =u   h ( x   1   , . . . ,x   d )=[ u   h   1   , . . . ,u   h   d ] T   εV   h   +u   0 ,  
         [0125]    such that  
                   ∫   Ω            ∑     i   ,   j   ,   k   ,     l   =   1       d            L   ijkl            ɛ   ij          (     u   h     )              ɛ   kl          (     w   h     )               x           =         ∫   Ω            f   ·     w   h               x         +       ∫     Γ   2              g   ·     w   h               s                     ∀       w   h     ∈     V   h                 ,           (   5   )                               
 
         [0126]    where the symbol “dot” denotes the standard vector multiplication.  
         [0127]    C.3 Problem-Oriented Criterion  
         [0128]    Users are often interested not only in the overall error e=u−u h , but also in its local behaviour, e.g., in a certain subdomain ∩⊂Ω.  
         [0129]    One way to get an information about the local behaviour of u−u h  is to measure the error in terms of specially selected problem-oriented criteria.  
         [0130]    One of the most typical criteria of such a type is presented by the integral  
                 ∫   Ω            Φ   ·     (     u   -     u   h       )               x         ,           (   6   )                               
 
         [0131]    where Φ=Φ(x 1 , . . . ,x d )=[φ 1 , . . . ,φ d ] T  is a selected vector-function such that Φ vanishes outside of the zone of interest (ω). C.4 Error Estimation  
         [0132]    see Section B4.  
         [0133]    C.5 Techniques for The Error Estimation in Terms of Problem-Oriented Criteria  
         [0134]    To present the techniques for estimation of criterion given by (6), we need first to describe two technical problems, that must be previously solved. They consist of finding an approximate solution of an auxiliary problem and making a certain post processing of solutions of this auxiliary problem and also of finite element approximation u h .  
         [0135]    C.6 Auxiliary Problem and its Finite Element Solution  
         [0136]    Let V r  be another finite-dimensional space constructed by means of a selected set of finite element trial functions on another standard finite element mesh T τ  over Ω. We notice that space V τ  is chosen so that its functions w τ  vanish on Γ, and also that T τ  need not to coincide with T h .  
         [0137]    Consider the auxiliary finite-dimensional problem as follows: Find a function  
           v   τ   =v   τ ( x   1   , . . . ,x   d )εV τ ,  
         [0138]    such that  
                 ∫   Ω            ∑     i   ,   j   ,   k   ,     l   =   1       d            L   ijkl            ɛ   ij          (     v   τ     )              ɛ   kl          (     w   τ     )               x           =       ∫   Ω            Φ   ·     w   τ               x                     ∀       w   τ     ∈       V   τ     .                     (   7   )                               
 
         [0139]    C.7 Tensor Averaging Procedures  
         [0140]    On T h , we define the tensor averaging transformation G h  mapping the tensor of the finite element approximation uh  
                 ɛ        (     u   h     )       =       [       1   2          (         ∂     u   h   i         ∂     x   j         +       ∂     u   h   j         ∂     x   i           )       ]       i   ,     j   =   1       d       ,           (   8   )                               
 
         [0141]    which is a constant d×d matrix over each element of the finite element mesh, into a d×d matrix-valued continuous piecewise affine function  
           G   h (ε( u   h ))=[ G   h   i,j (ε( u   h ))] i,j=1   d ,  (9)  
         [0142]    by setting each its nodal value as the mean value of ε(u h ) on all elements of the patch P(x o ) associated with corresponding node x o  in the mesh T h .  
         [0143]    More precisely, let x o  be one of nodes of finite element mesh T h , and let T 1 , . . . ,  
                     G   h     i   ,   j            (     ɛ        (     u   h     )       )            (     x   *     )       =       1     N     x   *                ∑       T   m     ∈     P        (     x   *     )                          1   2          (         ∂     u   h   i         ∂     x   j                       T   m            +       ∂     u   h   j         ∂     x   i                  T   m         )             ,     
        or           (     10      a     )                     G   h     i   ,   j            (     ɛ        (     u   h     )       )            (     x   *     )       =       1       measT   1     +   …   +     measT     N     x   *                    ∑       T   m     ∈     P        (     x   *     )                            measT   m     ·     1   2              (         ∂     u   h   i         ∂     x   j                       T   m            +       ∂     u   h   j         ∂     x   i                  T   m         )     .                   (     10      b     )                               
 
         [0144]    T N      x.    form the patch P(x o ) (cf. FIG. 2). Then, we define for i,j=1, . . . ,d  
         [0145]    Having G h   i,j (ε(u h )) defined at all nodes of the mesh T h , we uniquely define continuous piecewise affine function G h   i,j (ε(u h )) over the whole domain Ω. In this way, the matrix-valued continuous piecewise affine function G h (ε(u h )) in (9) is built. Similarly, on T τ , we define the tensor averaging transformation G τ  mapping the tensor of the finite element approximation v r  
                 ɛ        (     v   τ     )       =       [       1   2          (         ∂     v   τ   i         ∂     x   j         +       ∂     v   τ   j         ∂     x   i           )       ]       i   ,     j   =   1       d       ,           (   11   )                               
 
         [0146]    which is a d×d constant matrix over each element of the finite element mesh, into a d×d matrix-valued continuous piecewise affine function  
           G   τ (ε( v   τ ))=[ G   τ   i,j (ε( v   τ ))] i,j=1   d ,  (12)  
         [0147]    by setting each its nodal value as the mean value of ∇v τ  on all elements of the patch P(x o ) associated with corresponding node x o  in the mesh T τ .  
         [0148]    More precisely, let x o  be one of nodes of finite element mesh T τ , and let T 1 , . . . ,  
         [0149]    T N     xo    form the patch P(x o ). Then, we define for i, j=1, . . . ,d  
                   G   τ     i   ,   j            (     ɛ        (     v   τ     )       )            (     x   o     )       =         1     N     x   o                ∑       T   m     ∈     P        (     x   o     )                          1   2            (       ∂     v   τ   i         ∂     x   j                T   m             +         ∂     v   τ   j         ∂     x   i                     T   m       )                   or               (     13      a     )                     G   τ     i   ,   j            (     ɛ        (     v   τ     )       )            (     x   o     )       =       1       measT   1     +   ⋯   +     measT     N     x   o                    ∑       T   m     ∈     P        (     x   o     )                            measT   m     ·     1   2            (         ∂     v   τ   i         ∂     x   j                     T   m              +       ∂     v   τ   j         ∂     x   i                         T   m       )     .                           (     13      b     )                               
 
         [0150]    Having G r   i,j (ε(v τ ))(x o ) at all nodes x o , we uniquely define continuous piecewise affine function G τ   i,j (ε(v τ )) over the whole domain Ω. In this way, the vector-valued continuous piecewise affine function G τ (ε(v τ )) in (12) is built.  
         [0151]    C.8 Formula for Calculation of the Estimator  
         [0152]    Our method estimates the quantity (6) by the quantity E(u h ,v τ ) given by the following formula:  
           E ( u   h   ,v   τ ):= E   0 ( u   h ,v τ )+ E   1 ( u   h   ,v   τ ),  (14)  
         [0153]    where  
                   E   0          (       u   h     ,     v   τ       )       =         ∫   Ω            f   ·     v   f               x         +       ∫     Γ   2              g   ·     v   τ               s         -       ∫   Ω            ∑     i   ,   j   ,   k   ,     l   =   1       d            L   ijkl            ɛ   ij          (     u   h     )              ɛ   kl          (     v   τ     )               x               ,   and           (   15   )                   E   1          (       u   h     ,     v   τ       )       =       ∫   Ω            ∑     i   ,   j   ,   k   ,     l   =   1       d              L   ijkl          (         ɛ   ij          (     u   h     )       -       G   h     i   ,   j            (     ɛ        (     u   h     )       )         )            (         ɛ   kl          (     v   τ     )       -       G   τ     k   ,   l            (     ɛ        (     v   τ     )       )         )               x     .                   (   16   )                               
 
         [0154]    The functional E(u h ,v τ ) is directly computable once the approximations u h  and v τ  are found.  
         [0155]    C.9 Error Estimation Algorithm  
         [0156]    The algorithm is similar as that one presented in Section B9.  
       D—Advantages of the Technique Presented with Respect to Others Available in the Literature  
       [0157]    In the available scientific literature, several methods of computer verification of the accuracy of approximate solutions in terms of “problem-oriented” criteria has been presented. The main advantages of our method with respect to others can be summarized as follows:  
         [0158]    the estimator contains an extra term, E 0 , which is directly computable and contains major part of the error;  
         [0159]    the second term is effectively estimated by means of computationally cheap procedure of “gradient-averaging”;  
         [0160]    the estimator is valid for a wide spectrum of approximations including the case when the meshes for the original and auxiliary (adjoint) problems are different.  
       II Implementation in a System  
       [0161]    A Introduction  
         [0162]    Now we describe how the method can be implemented in hardware, software, firmware or combinations thereof, particularly in the form of engineering software. Actually, the nature of the proposed methodology allows a development of the program block “CHECKER”, FIG. 13, which can be easily incorporated in or added in some way to the most of the existing commercial software packages using finite element method as a main computational tool (e.g., Matlab®, ADINA®; ANSYS®, MATHCAD®). To illustrate that we shall use the PDE Toolbox of Matlab® package, especially designed for obtaining the numerical solutions for various classes of boundary-value and initial boundary-value problems.  
         [0163]    Thus, using the PDE Toolbox subroutines we have developed a new software, which can help the users of Matlab® not only to calculate numerical solutions for the boundary-value problems, but also effectively estimate their accuracy. Naturally all necessary routines could be integrated into one software package.  
         [0164]    B Software Product and its Usage  
         [0165]    [0165]FIG. 4 presents a layout of desktop of the software in a graphical user interface (GUI) and on a display. The desktop (control panel) contains six main blocks: “Define problem”—block  10 , “Mesh construction”—block  12 , “Reference solution”—block  14 , “Construction of adjoint problem”—block  16 , “Estimator”—block  18 , “Table”—block  20 . The software works in a standard PC having processing unit, program and data memory (RAM and hard disk), keyboard, mouse and display. The operating system is Microsoft® Windows, but the software can be implemented in other operating systems, too.  
         [0166]    Block  10  calls the PDE Toolbox main window (FIG. 5) for setting the problem: defining the solution domain, zone of interest, boundary conditions, coefficients, source function, which are exported later to our software, step A in FIG. 18.  
         [0167]    Block  12  is used for a construction of the first mesh using mesh specification parameters of PDE Toolbox of Matlab® such as maximum edge size, growth rate and PDE Toolbox of Matlab® grid generator. The block  12  calls FE-solver of PDE Toolbox for finding the first finite element approximation, step B.  
         [0168]    Block  14  is optional for finding the error (brute force method) if the user has time enough and computer resources for it.  
         [0169]    Block  16  is used to construct the second mesh and find the second finite element approximation for the adjoint problem, step C. It calls PDE Toolbox of Matlab® grid generator and its FE-solver for refining the second mesh in the zone of interest and also globally in the whole solution domain and computing the second approximation, step D.  
         [0170]    Blocks  12 ,  14 , and  16  have special fields, where the parameters of corresponding meshes such as number of triangles and number of nodes are displayed.  
         [0171]    Block  18  is designed to compute the estimator (E), step E. There is a possibility of choosing the type of gradient averaging and a parameter for calculating the integrals in the estimator (E), the number of triangles in the mesh used to calculate the estimator is shown in a special field. Steps 5-9 of paragraph B.9 form step E here.  
         [0172]    Block  20  is a table designed to show the results of calculations and values of different parameters used in the processes of calculations. If user considers an error being too large, he or she can recomputed with changed initial parameters (STEP F1 or F2).  
         [0173]    Referring now to FIG. 5, the solution domain and the zone of interest are set defining their geometry on the desktop of the PDE Toolbox. The coefficients and the source function are set in a known manner with the PDE Toolbox&#39;s facilities by using its panel buttons or drop menus.  
         [0174]    The solution domain presents here a thin plate, which is heated overall by a constant heat source (f=10 everywhere). The material of the plate is assumed to be homogeneous so that we can take the diagonal diffusion coefficients to be equal to 1, and the off-diagonal coefficients to be equal to 0. For simplicity, we assume that we have homogeneous Dirichlet boundary condition only, i.e., the temperature is assumed to be zero over the boundary of the plate. The temperature distribution in this plate can be described then by the following BVP  
             -         ∂   2        u       ∂     x   1   2           -         ∂   2        u       ∂     x   2   2           =     10                 in                 Ω       ,                         
 
         [0175]    The knowledge of the temperature, especially at some critical areas, is very important for making the right design solutions, which can save, e.g., material resources or help to improve the geometry of the plate. In the domains like we are dealing with now which have concavity in the boundary, such critical area is neighbourhood of concavity point, which dictated our choice of the zone of interest.  
         [0176]    After the setting of the problem, all the date of the problem are exported to our software, FIG. 6 shows PDE Toolbox&#39;s export window  26 . Parameters gd refers to “geometry data”, sf refers to “set formula”, and ns refers to “labels” 
         [0177]    [0177]FIG. 7 presents setting the boundary conditions using the Matlab® PDE Toolbox&#39;s means.  
         [0178]    [0178]FIG. 8 presents exporting the decomposed geometry and boundary conditions to CHECKER software using the export window  26 .  
         [0179]    The export functions used in FIGS. 6 and 8 are provided by Matlab® PDE Toolbox.  
         [0180]    After the previously defined data have been exported to the CHECKER software, the first finite element approximation (Primal Solution) is calculated using Block  12 . The result of calculations is shown in a new popup window  22 , which presents the 3D visualization of the approximate solution. The used mesh can be obtained (see FIG. 9) by rotating the 3D visualization.  
         [0181]    Block  14  by brute force method computes the reference solution on very fine mesh and tries to find the exact value of the problem-oriented criterion. This brute force method is not advisable for using in the engineering practice since it takes a large amount of computational time and memory. The brute force methods may require 100 times computer resources (processor speed, processing time and memory) as the new method herein described. The reference solution here is computed only for an illustration of the invention as it gives a more or less reliable reference result. The results of performance of Block  14  are shown in Table 20 in FIG. 12 (“exact error”). Ordinary user does not need Block  14 .  
         [0182]    Block  16  constructs the adjoint mesh and solves the adjoint problem (with PDE Toolbox). In the present implementation the function (ù) in the adjoint problem is set by default to 1 in the zone of interest and zero outside of it. Further implementation will make it possible to set the function (ù) arbitrarily, which can help to analyse critical areas more carefully. FIG. 10 shows the results of calculations of the adjoint problem in a new pop-up window  24  similarly to FIG. 9. In Block  16  we have also a possibility to construct another adjoint mesh condensed in the zone of interest and make new calculations, which is shown in 3D visualization in FIG. 11.  
         [0183]    After putting window  24  to the background we can make visible the table 20 with the results of all previous calculations (FIG. 12).  
         [0184]    All the blocks compute all the approximate solutions and the estimator using formulas as described above. Block  12  computes the first finite element approximation. Block  16  computes the second finite element approximation (for the adjoint/auxiliary problem). Block  18  computes the value of the estimator using the invented formula.  
         [0185]    After getting the result user is able to consider if the estimator shows that the first approximation isn&#39;t enough accurate, he/she can change initial parameters and recompute a new approximation and a new estimator.  
       III Model Example—One Dimensional Case  
       [0186]    The model depicts a heated rod, which is insulated but heated on its length. The rods ends are kept in zero temperature. The heating is defined by the source function. The designer gets advantage when he/she is able to compute temperature distribution more accurately and computationally cheaper in critical area.  
         [0187]    We begin with forming a simple problem: find the function u such that  
         − u {umlaut over ()}( x )= f ( x ) in Ω=(0,1),  (1)  
           u (0)= u (1)=0,  (2)  
         [0188]    where  
         f        (   x   )       =     {           2      α               exp        (       1     s   2       -     1   p       )       ·     [         -   2            (     x   -     x   0       )     2       +     4          p        (     x   -     x   0       )       2       +     p   2       ]     ·     1     p   4         ,                               if                 x     ∈     (     a   ,   b     )     ⋐     (     0   ,   1     )       ,               0   ,           otherwise   ,                                   
 
         [0189]    p=s 2 −(x−x 0 ) 2 , s=½(b−a), x 0 =½(a+b), and α is a positive real number.) Here a, b are coordinates of two sides of the rod. f(x) is a function that has a meaning of “heating function”. It shows t he amount of heat given to each elementary part of the rod. x is a coordinate changing along the rod, x 0  is a central point inside (a,b). s, p, and α: are positive parameters that define the shape of f(x). It should be noted that such a form of “f” is taken only for a particular example. u(x) has a meaning of the temperature of the rod at point x  
         [0190]    [0190]FIG. 14 presents the shape of the source function. The exact solution of problem (1)-(2) is  
         u        (   x   )       =     {             α   ·     exp        (       1     s   2       -     1   p       )         ,               if                 x     ∈     (     a   ,   b     )     ⋐     (     0   ,   1     )       ,               0   ,           otherwise   .                                   
 
         [0191]    The exact and approximate solutions of problem (1)-(2) (full line) and the zone of interest (dotted line) are marked in FIG. 15.  
         [0192]    As a quantity of interest, we take the following integral  
                   ∫   0   1            ϕ        (   x   )            (       u        (   x   )       -       u   h          (   x   )         )             x         ,     
        where          
            ϕ        (   x   )       =     {           1   ,             if                 x     ∈     (       z   1     ,     z   2       )                 0   ,           otherwise   .                       (   3   )                               
 
         [0193]    where  
         [0194]    In what follows, M denotes the number of elements in the mesh for the primal problem, N stands for the number of elements in the mesh used for the corresponding adjoint problem. In all the tests below, M=8, a=0, b=1.0, and α=1.0.  
         [0195]    1. Adjoint Problem is Solved on Uniformly Refined Meshes  
         [0196]    First, we present the results of numerical experiments where the uniform meshes are used, i.e., h=1/M and τ=1/N.  
       EXAMPLE 1  
       [0197]    In this series of numerical tests, (z 1 ,z 2 )=(0.375,0.625) (cf. FIG. 15).  
         [0198]    In Table I, we present the numerical results for problem (1)-(2) with M=8 and N=k·M, k=2,3,4,5,6,7,8.  
                                             TABLE I                           The results of performance of the estimator E for Example 1.            M   N   E 0     E 1     E   criterion   I eff                 8   16   0.00693883   0.00238174   0.00932058   0.00921705   1.01123239       8   24   0.00820668   0.00098981   0.00919649   0.00921705   0.99776967       8   32   0.00864915   0.00053743   0.00918658   0.00921705   0.99669427       8   40   0.00885372   0.00033653   0.00919025   0.00921705   0.99709278       8   48   0.00896479   0.00023026   0.00919505   0.00921705   0.99761349       8   56   0.00903174   0.00016737   0.00919911   0.00921705   0.99805336       8   64   0.00907518   0.00012710   0.00920229   0.00921705   0.99839820                  
 
         [0199]    The last column of the table presents the so-called “effectivity index”, which is often used to give a normalized measure of the quantity of error estimation. In our case, this index is computed as follows:  
               I   eff     =            E        (       u   h     ,     v   τ       )                     ∫   0   1            ϕ        (   x   )            (       u        (   x   )       -       u   h          (   x   )         )             x                        (   4   )                               
 
         [0200]    The dependence of I eff  on the number k is presented on FIG. 16. It is easy to see that the estimator E works quite well even for very modest values of N.  
         [0201]    2. Adjoint Problem is Solved on Condensed Meshes  
         [0202]    In this subsection, we present the results of numerical experiments, where the mesh used for solving the adjoint problem is condensed around the “zone of interest” (z 1 ,z 2 ). As before, h=1/M . Now, we are aimed to show that to have a good quality of the error estimation it is not necessarily required to also use globally refined meshes for the adjoint problem. We find that essentially the same result can be obtained on the mesh refined only in the “zone of interest” (i.e., we have less degrees of freedom and, thus, save computational costs). Thus, in the tests of Example 2, the mesh T τ  is only refined in a certain interval (q 1 ,q 2 ) that contains (z 1 ,z 2 ). To numerically qualify the respective meshes, we introduce positive numbers M i ,i=1,2,3, that correspond to the number of elements in the mesh for the primal problem in the intervals (0,q 1 ), (q 1 ,q 2 ), and (q 2 ,1.0), respectively. The numbers of elements in the mesh used for the adjoint problem in the intervals (0,q 1 ), (q 1 ,q 2 ), and (q 2 ,1.0) are defined in terms of two positive numbers k 1  and k 2  as follows—M 1 ·k 1 , M 2 ·k 2 , and M 3 ·k 1 , respectively. Thus, the density of the mesh used for the adjoint problem with respect to the mesh used for the primal problem inside the interval (q 1 ,q 2 ) is given by k 2 , and outside of the interval (q 1 ,q 2 )—by k 1 . Evidently, the total number of elements in the mesh used for the adjoint problem is given by N=M 1 ·k 1 +M 2 ·k 2 +M 3 ·k 1 , see FIG. 17, showing meshes T h  and T τ  used for problem (1)-(2).  
       EXAMPLE 2  
     We take (z 1 ,z 2 )=(0.500,0.750) and (q 1 ,q 2 )=(0.375,0.875). The results of the tests for various values of k i  are presented in Table II.  
       [0203]    [0203]                                                                                           TABLE II                           The results of performance of the estimator E for       Example 2.            M 1     M 2     M 3     k 1     k 2     M   N  E   criterion   I eff                      3   4   1   3   3   8   24 0.00495048   0.00500423   0.98925780       3   4   1   2   3   8   20 0.00495048   0.00500423   0.98925777       3   4   1   1   3   8   16 0.00495048   0.00500423   0.98925773       3   4   1   1   2   8   12 0.00493566   0.00500423   0.98629572       6   8   2   3   3   16   48 0.00122594   0.00122749   0.99873983       6   8   2   2   3   16   40 0.00122594   0.00122749   0.99873972       6   8   2   1   3   16   32 0.00122594   0.00122749   0.99873956       6   8   2   1   2   16   24 0.00123056   0.00122749   1.00249939                    
         [0204]    We see that the quality of the error estimation remains quite good even if  
           M   N     =   1.5     ,                         
 
         [0205]    what shows the effectivity of local mesh refinement in the “zone of interest”.  
         [0206]    Referring now to FIG. 12 to the results exposed in the table 20 we clearly observe good performance of the estimator for two dimensional problems. The “exact error” is based on the reference solution computed by brute force method used instead of the exact solution. The “exact error” means the relevant quantity of interest computed by the same formula (4) as presented paragraph B.3.  
         [0207]    In the above examples the models depict diffusion type processes, which arise often in various engineering applications. A typical diffusion process is heat propagation inside a body. So that physically our examples can be considered as problems on the distribution of a heat inside a body with given temperature at the boundary and a certain distribution (f) of the heat source. It is worth outlining that the algorithm suggested is applicable not only to such type of problems, but for all other physical models that can be described in terms of linear elliptic equations (e.g., linear elasticity, electrostatic problems, certain models in the theory of plates and shells). We refer to diffusion type problems as to a simple and transparent example.  
       REFERENCES  
       [0208]    [1]. Babu{haeck over (s)}ka, T. Strouboulis. The Finite Element Method and its Reliability, Oxford University Press Inc., New York, 2001.  
         [0209]    [2] R. Becker, R. Rannacher. A feed-back approach to error control in finite element methods: Basic approach and examples. East-West J. Numer. Math., 4, 237-264, 1996.  
         [0210]    [3] S. Korotov, P. Neittaanmäki, S. Repin. A posteriori error estimation of goal-oriented quantities by the superconvergence patch recovery. J. Numer. Math. 11, 33-59, 2003.  
         [0211]    [4] J. T. Oden, S. Prudhomme. Goal-oriented error estimation and adaptivity for the finite element method. Comput. Math. Appl., 41, 735-756, 2001.  
         [0212]    [5] R. Verfurth. A review of a posteriori error estimation and adaptive mesh-refinement techniques), Wiley-Teubner, 1996.