U’(Y)
= 
I
(2, g’(x, YKQ
f)(x) df?x
Ir,ii’k
Y)U~U%X) dT,
(45)
Note that the residuals of the computed coarsescale solution, that is 2&’ f and [&“], are the sources of error and the finescale Green’s function acts as the distributor of error. There seems to be agreement in the literature on a posteriori estimators that either one of, or both, the residuals are the sources of error. Where there seems to be considerable disagreement is in how the sources are distributed. From (45) we see that there is no universal solution to the question of what constitutes an appropriate distribution scheme. It is strongly dependent on the particular operator 2’ through the finescale Green’s function. This result may serve as a context for understanding differences of opinion which have occurred over procedures of a posteriori error estimation. REMARK. An advantage of the variational multiscale method is that it comes equipped with a finescale solution which may be viewed as an a posteriori estimate of the coarsescale solution error. DISCUSSION. Let us examine the behavior of the exact counterpart of (45) for different operators of interest. Assume that uh is piecewise linear in all cases, and that f= 0. First, consider the Laplace operator, 22 = A, b = a/&z. In this case, 2%” = 0, and the interface residual, [[b?l], is the entire source of error. Keeping in mind the highly local nature of the Green’s function for the Laplacian, a local distribution of [&“I would seem to be a reasonable approximation. The same could be said for linear elasticity, assuming there are no constraints, such as, for example, incompressibility, or unidirectional inextensibility. Next, consider the advectiondiffusion operator, .Z = a .V  KA, [[bi”‘] = [K dii”/dn~‘, where a is a given solenoidal velocity field, and K > 0, the diffusivity, is a positive constant. In the case of diffusion domination (i.e. advective effects are negligible), 2 = KA, and the situation is the same as for the Laplacian. On the other may be ignored. This time, the hand, when advection dominates, 2ipUh = a .VUh, and l[&“] = [K di?/dn] element residual, Z’ih, is the source of error. A local distribution scheme would seem less than optimal because the Green’s function propagates information along the integral curves of a (keep in mind that the Green’s function is for the adjoint operator, Z’*), with little amplitude decay. This means that there is an approximately constant trajectory of error corresponding to the residual error _c;PUh,in the element in question. Finally, consider the Helmholtz operator, 2 = A  k’, b = c~/&z, where k is the wave number. If k is real, we have propagating waves, whereas if k is imaginary, we have evanescent (decaying) waves. In the latter case, the Green’s function is highly localized; as Ikl +O the Green’s function approaches that for the Laplacian, as [kl+ cz the Green’s function approaches kP2S, a delta function. In the case of propagating waves, the Green’s function is oscillatory. In general, for Ikl large, the dominant source of error is the element residual, _Y? = k2Uh. As Ikl 0, the interface residual, [bU”] = [&“/&z], dominates.
4. Hierarchical
prefinement
and bubbles
Hierarchical prefinement plays an important role in clarifying the nature of the finescale Green’s function, g’, and provides a framework for its approximation. We begin by introducing some notations. Let
’ This follows
from the continuity
of advective
flux
T.J.R. Hughes et al. I Comput. Methods A&.
2” = c GALA (likewise A=,
w”)
(46)
where F’ is a finite element shape function corresponding nodal values; and let nn
Mrch. En,grg_ 166 (1998) 324
(likewise
with the primary nodes, A = 1,2, . . . ,n node,, and LA is the
associated
w’)
(47)
finite element shape function associated with the additional nodes, A = Ni is a hierarchical 1,2, , , Il;&,> and U; are the corresponding hierarchical degrees of freedom. For example, let Uh be expanded in piecewise linear basis functions and U’ in hierarchical cubits (see Fig. 9). Note, bubble functions are zero on element boundaries. Substituting (46) and (47) into (17)(20) and eliminating U; by static condensation results in where
B(w”,uh;g’)=L(Wh;g’)
VWhEClr
(48)
where B(Wh, Uh; g’) = a(Wh, 2)
+ (.L!?*G’, ti’( 25”))
L(w”; $) = (Wh, ,f) + (T*w”,
(49)
M’f)
(50)
and
(eY*i?)(y)g’(x,y)( 2&i”)(x)dR,
d0,
(51)
“nodr\
i’k
Y) = *:_, N;(Y)[(~"
(52)
'l,&,(x)
tih = standard linears (a) u’ = hierarchical cubits (o)
typical linear t
gal
,. typical cubic edge function
/
typical cubic
bubble function
Fig. 9. Hierarchical
cubits in two dimensions
T.J.R. Hughes et al. I Comput.MethodsAppl. Mech. Engrg. 166 (1998) 324
15
where
K" = [KA',]
(53)
K18 = a(NA, NL)
(54)
REMARKS (1) Recall, .2%Uhand 28’*Wh are Dirac distributions in the finite element case (cf. (51) and (39)). (2) Hierarchical prefinement generates an approximate finescale Green’s function, g’ = g’. in forms avoiding (3) For implementational purposes, it is more convenient to rewrite (49)(51) distributions. This can be accomplished by using the integrationbyparts formulas, viz.
Dirac
“no&a (Z*w”,
&i’( _‘SSh f))
= .T=,
a(Wh,N~>t(K”)‘l,,(a(N~, U”) WL, f)>
which amounts to the usual static condensation (4) A posteriori error estimation for the coarsescale (38) and (44)):
I =J‘ ==(I
u’(y) = 
n
i(x, YM'~~ Z(x, Y)Wih
R’
algorithm. solution, Uh, IS provided
by the finescale
(55)
solution
(see
f)(x) df4 f)(x) dfix  i,, ii’k ~>Ub~~ll(x)dT,
‘14
Rr g’(x, Y)@'~~
f>(x) dfX +
2’6, Y)@~~>(x>dT,
e=l
or, in analogy
with (55), by
U’(Y) =  A:=,
(5) (6)
(7)
(8)
(9)
(56)
>
N~(y>[(K”>‘l,,CaCN~,Uh> WL, fN
(57)
The quality of this estimator depends on the ability of {Ni}z$+s to approximate the fine scales, or equivalently, the quality of the approximation S’ = g’. Note that the finescale Green’s function only depends on the hierarchical basis. The exact finescale Green’s function corresponds to the limit p + 00. The finescale Green’s function is nonlocal, but it is computed from a basis of functions having compact support. For example, in the twodimensional case, the basis consists of bubbles, supported by individual elements, and edge functions, supported by pairs of elements sharing an edge. The threedimensional case is similar, but somewhat more complicated; the basis consists of bubbles, face and edge functions. In three dimensions, pairs of elements support face functions whereas the number of elements supporting edge functions depends on the topology of the mesh. In two dimensions, by virtue of the convergence of hierarchical prefinement, the exact finescale solution may be decomposed into a finite number of limit functionsone bubble for each element and one edge function for each pair of elements sharing an edge. In one dimension the situation is simpler in that only bubbles are required. The threedimensional case is more complex in that bubbles, face and edge functions are required. Polynomial bubbles are typically ineffective, but socalled residualfree bubbles [2] are equivalent to exactly calculating element Green’s functions. This approximation works exceptionally well in many important cases and will be discussed in more detail later on. If we consider onZy bubble functions in the refinement, coupling between different elements and all jump terms are eliminated. In this case,
(58)
T.J.R. Hughes et al. I
Cornput.Methods Appl. Mech. Gzgrg. 166 (1998) .?I4
where 0’ has replaced 0, and now the effect of u ’ is nonlocal only within and takes on zero approximate Green’s function g’, is defined elementwise boundaries.
5. Residualfree
each element. The values on element
bubbles
The concept of residualfree bubbles has been developed and explored in [l3,68,17,18]. The basic idea is to solve the finescale equation on individual elements with zero Dirichlet boundary conditions, viz. we wish to find u’ E Y’, such that Vu E z’, ZZ’_Yu’ = Zl’(%i
,f)
on R’
u’=O
on r’ >
e = I, 2,
, n2,,
(59)
Noting that u can be expressed in terms of the coarsescale basis having support in the element in question, realize that a finescale basis of residualfree bubbles can be constructed for each element, i.e.
we
where a = 1, 2, . , II,, is the local numbering of the primary nodes of element e. Thus, to each coarsescale basis function ??,,, we solve (60) for a corresponding residualfree bubble N:,. Consequently, the maximal dimension of the space of residualfree bubbles for element e is nen. It is typical, however, that the dimension is less than nen. Brezzi and Russo [3] have constructed residualfree bubbles for the homogeneous advectiondiffusion equation assuming the coarsescale basis consists of continuous, piecewise linears on triangles. For this case II = 3, but the dimension of the space of residualfree bubbles is only one. Let Be denote the residualfree b;bble basis solution of the following problem: Lf?B<, = 1
on R’
B, = 0
on r’ 1
(61)
Note, due to the fact the coarsescale space consist only of piecewise linears, combined with the fact that the finescale space satisfies zero Dirichlet boundary conditions, we can omit the projection operator, II’, present in the general case, namely (60) and solve (61) in the strong sense. However, in order to avoid potential linear dependencies, in general, we need to respect (60).
6. Element
Green’s functions
The idea of employing an element Green’s function was proposed in our initial work on the variational multiscale method [IO] (see also [ 131). In place of (23)(24) we solve the Green’s function problem for each element: WZ*g;,(X,
y) = n’s(x
 y)
g:cx> Y) = 0 Use of element
Green’s
g’(x,Y)=g:(x,y)
VXEr’
functions
onr’,
el,2
e=
1,2 ,...,
n,,
in place of the global Green’s e= 1,2 ,...,
scales vanish on element ,...,
(62)
1
vxx,yEfic,
The upshot is that the subgrid u’=O
V X E R’
function
amounts
to a locul upproximution,
n,, boundaries,
(63) i.e.
n,,
This means the subgrid scales are completely confined within element interiors. There is an intimate link between element Green’s functions and residualfree explored in [2], in which it was shown that, for the case governed by (61),
(64) bubbles.
This idea was first
T.J.R. Hughes
B,(Y) =
I
e
R
dk
et (11. I Comput.
Methods
Appl.
Mech.
Engrg.
166 (1998)
17
324
Y>df4
(65)
This result can be easily derived as follows:
I g:
da =
me
k:> zB,)w
=
= 4s:> BJtp = (Z*g:, BJnc, = (4 Be)p
= Be
(66)
Another way to derive (65) is to appeal to the general formula for the Green’s function basis, namely (52). Specialized to the present case, (52) becomes S:(X> Y)
= B&MB,,
in terms of a finescale
(67)
B,),,)‘B,(y)
Note that = (Be, ZB,),,
[email protected],> B&e
= (B,, l)ne
I
YE
By integrating
I
R 1’
(68)
Beda
R’
(67) and using (68), we have
g:,(x>Y) df4 = =
I
B,(x)dfU4B,, BJw'UY>
I>’
(69)
B,(Y)
REMARK. In general, the relationship result (65) is special to a residualfree
7. Stabilized Classical a$,
between a finescale basis and a Green’s function bubble governed by (61).
is given by (52). The
methods stabilized
methods
U”) + (LWh, 7( Z?
where L is typically L = +2? L = +Tudv L = 2?*
are generalized f))n,
a differential
Galerkin/
methods
of the form
= (Wh, f)
operator,
leastsquares
Galerkin
(70)
such as
(GLS)
(71)
SUPG
(72)
Adjoint
(73)
and r is typically an algebraic operator. SUPG is a method defined for advectivediffusive decomposable into advective (d%j,,) and diffusive (.Z?&) parts. 7.1. Relationship
of stabilized
methods
with subgridscale
operators,
i.e. ones
models
It was shown in [lo] that a stabilized method of adjoint type is an approximate subgridscale model in which the algebraic operator 7 approximates the integral operator M’ based on element Green’s functions,
18
T.J.R. Hughes
et 01. I Comput.
Methods
Appl.
Mech.
Engrg.
166 (1998)
324
T=il;JIIM’
(74)
Equivalently, T. S(y  x) = g:tx, y) = g’k
L’)
II
25”
(75)
Viz.
R’
(2’*?)(
y,g:,(x, y)(
f)(x)
dJ& dfJ,(2,.
0’
(c%‘*W”)(y)~. 6( y  x)(.25?.f)(x> dR, df2,
=
(T*i”)(x)~ . (cYi? ,f)(x)
Lx
dR,
(76)
7.2. Formula for r based on the element Green’s function The approximation 2. &Y

x) =
suggests defining
g:cx> Y)
(77)
r by (78)
(79)
REMARKS (1) The element mean value of the Green’s function provides the simplest definition of r. (2) This formula is adequate for loworder methods (hadaptivity). For higherorder methods (padaptivity), accounting for variation of r over an element may be required. In this case, we may assume for example that T = 7(x, y) is a polynomial of sufficiently high degree. Given an element Green’s function g: we can in principle always calculate an equivalent function r. Consequently, there is a generalized stabilized method, i.e. a method of the form, a(Wh, 2) 
cl<
*’ I *e
2Y*Wh( )‘)T(X, y)( 2%” f)(x)
df2, d.ny = (Wh, f‘>
equivalent to the element Green’s function method. The generalized determining r such that the following equivalence condition is satisfied
JI 0“ RC
L?*W”( y)~(x, y)( 3;”
f)(x)
stabilized
(80) method
involves
dac df$
Thus, we have a full equivalence as indicated procedure will be presented next.
in Fig. 10. Examples
of the calculation
of 7 by way of this
EXAMPLES. We consider two onedimensional examples, the advectiondiffusion equation and the Helmholtz equation. In both cases we assume the use of standard piecewise linears for the coarsescale basis. Consequently, determination of the element Green’s function may be performed using the strongform counterpart of (62), namely
T.J.R. Hughes et al. I Comput. Methods Appl.
Mech. Engrg. 166 (1998) 324
19
element Green’s functions
residual free bubbles <_
stabilized methods Fig. 10. Stabilized methods can be constructed which are equivalent to methods based on element Green’s functions, which in turn are equivalent to methods developed from the residualfree bubbles concept. Historically, stabilized methods preceded the residualfree bubbles and element Green’s function approaches, and there are certain stabilized methods (e.g. SUPG and GLS) which have been established for some time that are not, strictly speaking, equivalent to these concepts. Nevertheless, they have been justified independently by mathematical analysis and numerical testing.
.L?*g:(x, y) = S(x  y)
vx
g:cG Y) = 0
VXEI’
REMARK. for any
In one dimension,
E
ne
e = 1,2,.
. , n,,
(82)
use of the element Green’s function
g: results in II =U + U’ being pointwise
exact
(83) and u being an endnode
interpolant
7.2.1. Advectiondiffusion Let
equation
of u. This result holds for all problems. ([IO])
(84) where .Z& = a i
(85)
$
Tdlff = K and a and
K
(86)
are assumed
Zu =f
to be positive
constants.
in 0 = [0, L]
u=O
Consider
the homogeneous
Dirichlet problem (87)
on r={O,L}
(88)
Note that _Y* = X$,
+ Z&r + Ziiff
= %dv The solution
(89) (90)
of (82) is given by
g:(A Y) =
1
C,(y)(l
 e~2ar’h)
C,(y)(ep2”‘“‘h’ em’“)
XSY x>y
(91)
where le
C,(Y)
=
C,(Y) = e
2u(
1 (y/h))
a( 1  ee2”) Za(ylh) 1
(93)
a(1  e*“)
ah (Y = 2~
(element
(92)
P&let number)
(94)
20
T.J.R. Hughes et (11. I Camput.
:d)
Fig.
0
=
5
I I. Element Green’s function for the onedimensional
and h = meas P&let numbers. By virtue of the each element, the use of the element
is the element
Methods
fe)
(2
Appl.
=
Mwh.
10
,i:
advectiondiffusion
length. This element
Eqr,q.
166 (1998)
u =
equation
Green’s function
324
i”O
as a function
of element PCCKL numoer t(u)
is shown in Fig. 11 for various element
fact that the coarse scales are piecewise linear and 2;’ and Z*w” are therefore constants on simple formula (79) suffices to define a constant S satisfying (Xl), that is, one equivalent to Green’s function, viz.
(9S) REMARKS ( 1) This is a wellknown formula from the theory of stabilized methods. (2) This r results in a nodully exact stabilized method for piecewise linear GA’s and c0nstant.f. element lengths need not be uniform. (3) The approximate fine scale solution is given by
Remarkably,
For the advectiondominated case we can neglect the jump term and make use of the approximation L!?Zh =a .VUh. Further, employing the approximation (75), we have u’ E [email protected] .VU” f‘) Assuming
u’ =
(97)
r has the form h /( 21ul) in the advectiondominated  +a 21al
.VUh f)
This may be used as a simple,
local error estimator
7.2.2. Helmholtz equation (OberaiPinsky /16/) The setup is identical to the previous example. Helmholtz equation in which
_Cf,!$2
case, (97) becomes
for the advectiondominated
We consider
the Dirichlet
case.
problem
(87)(S),
for the
(99)
is the Helmholtz operator and k is the wave number. We assume k E R, corresponding to the case of propagating waves. Note Z’* = _Y. The Green’s function for an element of length h is given by
T.J.R. Hughes et al. I Comput. Methods Appl. Mech. Engrg. 166 (1998) 324
(b) kh = VT
(aj kh = 2n/lO
(d)
(c) kh = 217
Fig. 12. Element Green’s
function
(kg:) for the onedimensional
Helmholtz
kh = 4r
equation
as a function
of element wave number (kh).
sin(kh( 1 + 5) / 2) sin(kh( 1  7) / 2) 3 5
k sin(kh)
g’(s’rl)=
sin(kh(1  5)/2)sin(kh(l k sin(kh)
where 5 and v are normalized,
biunit
(100)
+r])/2) 3 527 coordinates,
1 a&,77c+1
(101)
The element Green’s function is depicted in Fig. 12 for various element wave numbers. This time .Zih and Z*Wh vary linearly over each element and satisfaction of the equivalence entails a nonconstant 7, as follows: 7( & rl) = 70” + r1,srl=
g:t &rl)
condition,
(81)
(102)
where
(103)
REMARKS (1) Oberai and Pinsky [16] have also derived an rectangle and an equivalent r. Similar results (2) As in the previous example, this r results in and constant J for an arbitrary nonuniform
8. Beyond residualfree
element Green’s function in two dimensions for a bilinear for three dimensions are also easily derived. a nodully exact stabilized method for piecewise linear &‘s mesh.
bubbles
The residualfree bubbles concept is established by way of its equivalence with mathematically proven stabilized methods. To make further improvement, we need to acknowledge the presence of functions other than bubbles. As discussed in Section 4, bubbles and edge functions can be constructed in two dimensions in such a way as to decompose the exact finescale solution in terms of functions with local support. To determine the
22
T.J.R. Hughes et al. I Comput. Methods Appl. Mech. Engrg.
166 (1998) 324
exact basis of bubbles and edge functions is problematic due to the fact that they do not have disjoint support. Consequently, uncoupled, local evaluation is not feasible and the coupled analytical problem appears virtually impossible to solve. However, it would seem that improvement, beyond that achieved with bubbles alone, could be attained by adding some kind of edge functions to the formulation to reduce the distributional component of error concentrated on element interfaces. This leads us to consider approximate problems governing residualfree functions in which simplifications have been introduced ab initio to avoid coupling. As an example, consider the setup for two contiguous linear triangular elements illustrated in Fig. 13. The idea is to define residualfree edge functions on pairs of elements ignoring the coupling introduced by edge functions associated with overlapping pairs of elements. Intuitively, functions defined in this way would appear beneficial but, of course, not capable of removing the distributional part of error entirely. For each pair of contiguous elements, the following problem is to be solved: Find Ni E w“(w), A = 1,2, 3,4, such that V u’ E Y(w), (u’, ZN:,)W = (u’,
[email protected], + B,
+ B;)
= (u’,
z(??, + B,
+ B,t) f),,
where Ni is a residualfree and
edge function,
Y’(w) = {u EH’(w)lu
= 0
f),,,
B,
 (u’, [b(i,
+ B, + B,‘)]),,
and B,f are bubbles
with support in R
on r}
(104) and R’,
respectively,
(105)
(Due to possible linear dependencies, the number of independent bubbles and edge functions may be reduced.) To make things a bit more specific, we can select bubbles to remove residuals on element interiors, as follows: (u’, c!ZB; )fI
= (u’,
(u’, Z&T+) f2+ = (u’,
L?&
f)n,
Vu’ E
d;p%<,+,f‘),),
V’(Li) ,
Vu’ E ‘V”‘(f2+),
a = 1,2, 3
(106)
a+ = 1, 2, 3
(107)
where
V’(C)
= {u E H’(fT)Ju
= 0
on r}
(108)
V’(0’)
= {u E H’(fl+)lu
=0
on r’}
(109)
Then, the edge functions (u’, ZN:,),
= (u’,
need to satisfy [[b(%, + B,
+ B,‘)])ys
Vu’ E ‘V(w),
Note that, in this case, ZNL = 0 on o’, but NI, is nontrivial Further simplifications may be envisioned. For example, the bubbles in (104), and define the edge functions by
A = 1, 2, 3,4
due to the jump term on y’ which acts as a source. one could ignore the distributional contribution of
JGJ3 g:::i a = 1
r=ifKl
fJ=nun+ w
=gs
y=rnT‘+ 7
=
(r u
(110)
r+)\y
Fig. 13. Typical patch of two contiguous linear triangles.
T.J.R. Hughes et al. I Comput. Methods Appl. Mech. Engrg.
(u’, Bv:,)w = (II’, Various other possibilities future work.
[bzJl),,
VU’ E v’(w),
suggest themselves.
166 (1998) 324
(111)
A = 1,2,3,4
We hope to pursue
23
this topic by way of specific examples
REMARK. These ideas are related to recent work of Franca and Russo [5] in which the residualfree is satisfied on macroelements.
in
condition
9. Summary We have given a general variational multiscale formulation of an abstract Dirichlet problem. We considered the case in which finite elements are used to represent the coarse scales and we derived an exact equation governing the coarse scales. This equation amounts to the Gale&in formulation plus an additional term which depends on the distributional form of the residual (i.e. element interior and interface jump terms) and a tinescale Green’s function, A representation is also derived for the finescale solution, which amounts to the error in the coarse scales. The results serve as paradigms for subgridscale models and a posteriori error estimators. To understand the nature of the finescale Green’s function, we considered hierarchical prefinement and bubbles. An explicit formula for the finescale Green’s function is obtained in terms of the hierarchical (i.e. fine scale) basis. This formula suggests that the finescale Green’s function can be represented in terms of a basis of functions having local support. Subsequent discussion dealt with developing practical approximations. We reviewed the concepts of residualfree bubbles, element Green’s functions, and stabilized methods, and summarized the relationships between them. We also presented some initiatory ideas on developing betterquality, approximate, residualfree bases by including edge functions in addition to bubbles in two dimensions, etc. We believe the variational multiscale method has considerable potential for important classes of problems in science and engineering. The major obstacles to be faced are explicit, analytic development of better quality, residualfree bases for problems of practical interest. Various ideas are being pursued at this time to address these issues.
Acknowledgements Brazilians and Italians are almost as good at soccer as they are at mathematics. We acknowledge the inspirations provided by our collaborators Franc0 Brezzi, Leo Franca and Alessandro Russo. This research was supported by the Office of Naval Research under Grant N000149610109.
References [I] C. Baiocchi, F. Breazi and L.P. Franca, Virtual bubbles and the Galerkin/leastsquares method, Comput. Methods in Appt. Mech. and Engrg. 10.5: (1993) 125142. [2] F. Brezzi, L.P. Franca, T.J.R. Hughes and A. Russo, h = 5 g, Comput. Methods in Appl. Mech. and Engrg. 145 (1997) 329339. [3] F. Brezzi and A. Russo, Choosing bubbles for advectiondiffusion problems, Math. Models Methods Appl. Sci. (1994) 571587. [4] J. Douglas Jr. and J.P. Wang, An absolutely stabilized finite element method for the Stokes problem, Math. Comput. 52 (1989) 495508. [5] L.P. Franca and A. Russo, Approximation of the Stokes problem by residualfree macro bubbles, EastWest J. Numer. Anal. 4 (1996) 265278. [6] L.P. Franca and A. Russo, Deriving upwinding, mass lumping and selective reduced integration by residualfree bubbles, Appl. Math. Lett. 9 (1996) 8388. [7] L.P. Franca and A. Russo, Mass lumping emanating from residualfree bubbles, Comput. Methods Appl. Mech. Engrg. 142 (1997) 353360. [8] L.P. Franca and A. Russo, Unlocking with residualfree bubbles, Comput. Methods in Appl. Mech. Engrg. 142 (1997) 361364.
24
T.J.R. Hughes
[9] L.P. Franca,
[ 101
[I I]
[ 121 [I31 [I41 [IS] [I61
[ 171 1181
et al.
I Comput.
Methods
Appl.
Mrch.
Engrg.
166 (1998)
324
T.J.R. Hughes and R. Stenberg, Stabilized finite element methods for the Stokes problem, In: M.D. Gunzburger and R.A. Nicolaides, eds., Incompressible Computattonal fluid dynamic\ (Cambridge University Press, Cambridge, 1993) 87107. T.J.R. Hughes, Multiscale phonomena: Green’s functions, the DirichlettoNeumann formulation, subgrid scale models, bubbles and the origins of stabilized method\, Comput. Methods Appl. Mech. Engrg. 127 (1995) 3X7401. T.J.R. Hughes and L.P. Franca, A new tinite element formulation for computational fluid dynamics. VII. The Stokes problem with various wellposed boundary conditions: Symmetric formulations that converge for all velocity/pressure spaces, Comput. Methods Appl. Mech. Engrg. 65 (1987) 8596. T.J.R. Hughes and G.M. Hulbert, Spacetime finite element methods for elastodynamics: Formulation5 and error estimates. Comput. Methods Appl. Mech. Engrg. 66 (1988) 339363. T.J.R. Hughes and J. Stewart, A spacetime formulation for multiscale phenomena, J. Comput. Appl. Math. 74 (1996) 217229. G.M. Hulbert and T.J.R. Hughes, Spacetime finite element methods for secondorder hyperbolic equations. Comput. Methods Appl. Mech. Engrg. 84 (1990) 3277348. K. Jansen, C. Whiting, S. Collis and F. Shakib, A better consistency for loworder stabilized tiniteelement methods, Preprint, 1997. A.A. Oberai and P.M. Pinaky, A multtscale finite element method for the Helmholti equation, Comput. Methods Appl. Mech. Engrg. I54 (1998) 2X1297. A. Russo, Bubble stabilization of tinite element methods for the IineariLed incompressible NavierStokes equations. Comput. Methods Appl. Mech. Engrg. 132 (1996) 335343. A. Russo, A posteriors error esttmators via bubble functions, Math. Models Methods Appl. Sci., 6 (1996) 3341.