2003 Summer Project

An hp-adaptive Finite Element (FE) Method for Solving Electromagnetic (EM) Problems. Part I: Petroleum Engineering applications. David Pardo Supervisor: A. Bespalov

Baker Hughes 2001 Rankin Rd. Houston, Tx 77073-5114 Abstract A fully automatic hp-adaptive strategy has been recently developed at ICES (Institute for Computational Engineering and Sciences) at The University of Texas at Austin. The method produces automatically a sequence of optimal hp-meshes that deliver exponential convergence rates for a large class of problems, including EM problems with some types of singularities. On this report, we present the adaptive strategy in context of Maxwell’s equations, and we study the suitability and utility of the method for Petroleum Engineering applications.

Acknowledgment To Dr. Bespalov, Dr. Demkowicz, Dr. Tabarovski and Dr. Torres-Verdin. Their help has been essential in order to develop the research presented in this report. 1

Contents 1 Introduction

4

2 Maxwell’s equations

5

2.1

Time harmonic Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . .

5

2.2

Electrostatics and Magnetostatics . . . . . . . . . . . . . . . . . . . . . . . .

6

2.3

Axisymmetric Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . .

6

2.4

Maxwell’s equations with Harmonic Dependence upon the Transverse Spatial Variable z . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

Boundary Conditions (BC’s) . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.5

3 Variational Formulation

10

3.1

Variational Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

3.2

Regularized Variational Formulation . . . . . . . . . . . . . . . . . . . . . .

11

3.3

Lagrange Multipliers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

4 A Fully Automatic hp-adaptive Finite Element (FE) Method

14

4.1

hp-Finite Elements (FE). . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

4.2

hp-edge Finite Elements. The de Rham diagram. . . . . . . . . . . . . . . . .

14

4.3

Survey of hp- and hp-edge FE codes. . . . . . . . . . . . . . . . . . . . . . .

15

4.4

2Dhp90 (3Dhp90): An anisotropic 2D (3D) fully automatic hp-adaptive package 16

4.5

The fully automatic hp-adaptive strategy

. . . . . . . . . . . . . . . . . . .

16

4.6

Goal-oriented hp-adaptivity . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

4.7

The need for a two grid solver . . . . . . . . . . . . . . . . . . . . . . . . . .

17

5 Analysis of Edge Singularities Arising in Electrostatic Computations

18

5.1

Model Problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

5.2

Numerical Results Corresponding to Model Problem 1: 270 degrees . . . . .

21

5.3

Numerical Results Corresponding to Model Problem 2: 358 degrees . . . . .

22

6 Limitations of the fully automatic hp-adaptive strategy 6.1

An axisymmetric battery in free space . . . . . . . . . . . . . . . . . . . . .

2

25 28

6.2

hp-edge Finite Element approximations using the fully automatic hp-adaptive strategy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

7 Summary and Conclusions

34

8 Future Work

35

9 Appendix: de Rham diagram

37

References

39

3

1

Introduction

Nowadays, there exist a variety of electromagnetic (EM) applications, ranging from highspeed communications and computing up to medical applications (see [30] for a detailed list). All these applications touch our daily lives and, therefore, they deserve a detailed study. A particular family of EM applications is of interest to us for this research: Petroleum Engineering applications. More precisely, the modeling of Logging While Drilling (LWD) EM measuring devices. In recent years, more advanced LWD tools have been designed. The use of LWD tools has become popular, since they allow for the correction of drilling trajectories in real time. These tools incorporate a number of measuring devices such as arrays of acoustic antennas, electromagnetic coils, etc. But due to the increasing number of measuring devices in a LWD tool, they are also geometrically more complex. Thus, the analysis of the interaction of EM fields with metallic structure of the LWD tool and surrounding areas cannot be solved by using analytical techniques only, and numerical methods are needed. These numerical methods should model details of the tool geometry, and provide an accurate answer in the area of interest (receiving coils). For details, see [20]. In this context, it is desirable a numerical technique that provides low discretization errors and, simultaneously, solves the discretized problem without prohibitive computational cost. Our preferred choice is a fully automatic hp-adaptive strategy developed by Dr. Demkowicz’s team at The University of Texas at Austin. Throughout this report, we analyze advantages and limitations of using this adaptive algorithm when solving Petroleum Engineering EM problems. In Section 2, we present Maxwell’s equations, which govern the EM phenomena. The corresponding variational formulation is derived in Section 3. In Section 4, the fully automatic hp-adaptive algorithm described in [14] is outlined. This adaptive strategy will be used in order to solve 2D (two dimensional) EM problems with edge singularities (see Section 5). In Section 6, limitations of the fully automatic adaptive strategy are illustrated with a real life Petroleum Engineering application. Section 7 is devoted toward conclusions whie guidelines of future work are suggested in Section 8.

4

2

Maxwell’s equations

All electromagnetic applications of interest for us are governed by Maxwell’s equations on a given domain Ω with prescribed material data and boundary conditions. Maxwell’s equations in the transient form are a set of two vector valued first-order partial differential equations (PDE’s) formulated in terms of the electric and magnetic fields: 1. ∇ × E = − 2. ∇ × H =

∂µH ∂t

∂²E + σE + Jimp ∂t

where E denotes the electric field, H the magnetic field, ² the permittivity of the medium, µ the permeability of the medium, σ the conductivity of the medium, Jimp the impressed electric current, and t the time.

2.1

Time harmonic Maxwell’s equations

In our applications, sources (antennas) are time harmonic with given radial frequency ω. Fourier transformation provides an equivalent system of equations, called time harmonic Maxwell’s equations, given by: ˜ = −jµω H ˜ 1. ∇ × E ˜ = jω²E ˜ + σE ˜ +J ˜imp 2. ∇ × H where ˜ indicates Fourier transform (phasors). For simplicity in notation, we will avoid the use of the symbol ˜ in the remaining of this report. Under certain regularity assumptions [16], this system of two first-order PDE’s can be reduced to two decoupled second-order PDE’s: one for the electric field, and another for the magnetic field, Ã

!

1 ∇ × E − (ω 2 ² − jωσ)E = −jωJ imp ∇× µ Ã

!

1 1 ∇× ∇ × H + jωµH = ∇ × J imp jω² + σ jω² + σ

5

(2.1)

(2.2)

Existence and Uniqueness of solution. Equations (2.1) and (2.2) have a unique solution if and only if ω does not vanish (transient case). In the steady state case (ω = 0), solution is unique up to a gradient. That is, uniqueness of solution will be guaranteed if and only if an extra condition on gradients is imposed.

2.2

Electrostatics and Magnetostatics

In the steady state case, Maxwell’s equations are reduced to the following system, (

∇×E =0 ∇ × H = σE + J imp

(2.3)

As mentioned above, solution of this system of equations is not unique, and a condition on gradients should be imposed. Electrostatics. If σ = 0, the use of Gauss law provides the following system: (

∇×E =0 ∇ · (²E) = ρ

(2.4)

where ρ is the distribution of free charges. In the case of a conductor (σ > 0), the free charges move, and we cannot prescribe them. Nevertheless, electrostatic equations can be derived by combination of the continuity equation and Faraday’s law: (

∇×E =0 ∇ · (σE) = ∇ · J imp

(2.5)

Magnetostatics. Once electric field is known, magnetostatic equations can be derived using system (2.3) together with Ampere and Gauss laws: (

2.3

∇ × H = σE + J imp ∇ × (µH) = 0

(2.6)

Axisymmetric Maxwell’s equations

In the area of Petroleum Engineering, it is commonly assumed that computational domain has an axial symmetry, i.e., both geological formation as well as logging tools are axisymmetric. Thus, extra attention should be given to this type of problems. 6

If the problem is axisymmetric, one can prove using symmetry considerations that both electric and magnetic fields are independent of the azimuthal variable φ. Thus, using cylindrical coordinates (ρ, φ, z), and assuming no dependence of (electric or magnetic) vector field A upon φ, we obtain: ∇ × A = −ˆ ρ

1 ∂(ρAφ ) ∂Aφ ˆ ∂Aρ ∂Az + φ( − ) + zˆ ∂z ∂z ∂ρ ρ ∂ρ

And Maxwell’s equations can be reduced to a system of six scalar valued first order partial differential equations (PDE’s):     (1)          (2)             (3)     (4)          (5)            (6)

∂Eφ = −jωµHρ ∂z ∂Eρ ∂Ez − = −jωµHφ ∂z ∂ρ 1 ∂(ρEφ ) = −jωµHz ρ ∂ρ ∂Hφ − = (jω² + σ)Eρ + Jρimp ∂z ∂Hρ ∂Hz − = (jω² + σ)Eφ + Jφimp ∂z ∂ρ 1 ∂(ρHφ ) = (jω² + σ)Ez + Jzimp ρ ∂ρ −

(2.7)

This system of six equations decouples into the following two systems:     (1)       

(3)

and

        (5)       (2)      

(4)

         (6)

∂Eφ = −jωµHρ ∂z 1 ∂(ρEφ ) = −jωµHz ρ ∂ρ ∂Hρ ∂Hz − = (jω² + σ)Eφ + Jφimp ∂z ∂ρ

(2.8)

∂Eρ ∂Ez − = −jωµHφ ∂z ∂ρ ∂Hφ = (jω² + σ)Eρ + Jρimp − ∂z 1 ∂(ρHφ ) = (jω² + σ)Ez + Jzimp ρ ∂ρ

(2.9)

−

7

2.4

Maxwell’s equations with Harmonic Dependence upon the Transverse Spatial Variable z

Now, we assume that electric and magnetic fields are harmonic with respect spatial variable z. Then, original 3D Maxwell’s equations can be reduced to a 2D problem for each z-frequency, leading to the so called 2.5D problems. By doing so, CPU time and memory required to solve the problem will decrease considerably. In this subsection, we present Maxwell’s equations assuming harmonic dependence upon variable z. Under this assumption, we obtain for a vector (electric or magnetic) field that: ∇ × A = xˆ(

∂Az ∂Ay ∂Ax ∂Az − jζAy ) + yˆ(jζAx − ) + zˆ( − ) ∂y ∂x ∂x ∂y

where ζ is the z-frequency And Maxwell’s equations can be reduced to a system of six scalar valued first order PDE’s:     (1)           (2)           (3)      (4)           (5)           (6)  

∂Ez − jζEy = −jωµHx ∂y ∂Ez = −jωµHy jζEx − ∂x ∂Ey ∂Ex − = −jωµHz ∂x ∂y ∂Hz − jζHy = (jω² + σ)Ex + Jximp ∂y ∂Hz = (jω² + σ)Ey + Jyimp jζHx − ∂x ∂Hy ∂Hx − = (jω² + σ)Ez + Jzimp ∂x ∂y

(2.10)

Remark: If ζ = 0, then (2.10) decouples into two systems of three scalar valued first order PDE’s, as in the case of axisymmetric problems.

2.5

Boundary Conditions (BC’s)

A large variety of BC’s can be imposed to EM problems. In this subsection, we briefly discuss the most commons ones, corresponding to Dirichlet, Neumann and Cauchy BC’s: • Dirichlet (essential) BC: 8

Homogeneous: It is used to model a perfect electric conductor (PEC). Since E in a PEC vanishes and it must be continuous across material interfaces, we obtain: n×E =0 Non-homogeneous: For scattering problems, E = E inc + E sca . Then, n × E sca = −n × E inc • Neumann (natural) BC: It may be used to model an antenna by prescribing an impressed surface current Jimp on the tangential component of the magnetic field. In terms of electric field E, it is given by: n × (∇ × E) = −jωµJimp S • Cauchy BC: Also known as impedance BC, it may be used to model E reflection from a material with large (but finite) conductivity: n × (∇ × E) + jωµβ(n × E) = −jωµJimp S

9

3

Variational Formulation

In this section, we first derive the variational formulation corresponding to Maxwell’s equations. Then, extra terms are added into the variational formulation in order to control gradients explicitly. As we will show, these extra terms are introduced through a Lagrange multiplier (also called scalar potential), which becomes crucial for steady state problems, as well as for problems with harmonic dependence upon a spatial variable. For simplicity, in the remainder of this section we consider only the reduced wave equation corresponding to the electric field. Notice that analogous formulations can be derived for the magnetic field reduced wave equation. Also, computational domain Ω is assumed to be simply connected, avoiding the technical issues associated with cohomology spaces, see e.g. [5].

3.1

Variational Formulation

In this subsection, we outline the derivation of the variational formulation corresponding to Maxwell’s equations. For details on the subject, see [8]. First, we multiply eq. (2.1) by the complex conjugate of test function F. Then, we integrate over computational domain Ω. Finally, integration by parts yields to the following formula: Z Z Z 1 ¯ )}dx + ¯ t dS = −jω J imp F ¯ dx { (∇ × E)(∇ × F n × (∇ × E)F Ω µ ∂Ω Ω At this point, we must introduce BC’s into the formulation. For the sake of argument, we assume ∂Ω = Γ1 ∪ Γ2 ∪ Γ3 , with a homogeneous Dirichlet BC on Γ1 , a Neumann BC on Γ2 and a Cauchy BC on Γ3 . Then, our variational formulation becomes:   Find E ∈ HD (curl; Ω) such that      Z Z   1        

Ω

µ

¯ )dx − (∇ × E) · (∇ × F −jω

½Z

Ω

Ω

¯ dx − jω (ω ² − jωσ)E · F 2

¯ dx + J imp · F

Z

Γ2 ∪Γ3

¯ dS J imp ·F S

Z

¾

Γ3

¯ dS = β(n × E)F for all F ∈ HD (curl; Ω) . (3.11)

In the above, HD (curl; Ω) is the Hilbert space of admissible solutions, H D (curl; Ω) := {E ∈ L2 (Ω) : ∇ × E ∈ L2 (Ω), n × E = 0 on Γ1 } .

(3.12)

with inner product defined by: (u, v)HD (curl;Ω) := (u, v)L2 (Ω) + (∇ × u, ∇ × v)L2 (Ω) 10

(3.13)

The original and variational formulations are equivalent to each other.

3.2

Regularized Variational Formulation

It is well known that the standard variational formulation is not uniformly stable with respect to wave number k 2 = µ(ω 2 ² − jωσ). As k → 0, we loose the control over gradients. This corresponds to the fact that, in the limiting case k = 0, the problem is ill-posed as the gradient component remains undetermined. A remedy to this problem is to enforce the 1 (Ω) continuity equation explicitly at the expense of introducing a Lagrange multiplier p ∈ H D of scalar potentials 1 HD (Ω) := {q ∈ H 1 (Ω) : q = 0 on Γ1 }, (3.14) 1 (Ω), to learn that solution E to Then, we employ a special test function F = ∇q, q ∈ HD (3.11) must automatically satisfy the weak form of the continuity equation,

−

Z

2

Ω

(ω ² − jωσ)E · ∇¯ q dx = −jω

½Z

J

imp

Ω

· ∇¯ q dx +

Z

Γ2

J imp S

· ∇¯ q dS

¾

.

(3.15)

We also recall the Helmholtz decomposition: E = ∇φ + E 0 ,

1 1 where φ ∈ HD (Ω) and (E 0 , ∇q)L2 (Ω) = 0 ∀q ∈ HD (Ω)

(3.16)

Then, the so called regularized variational formulation looks as follows:  1  Find E ∈ HD (curl; Ω), p ∈ HD (Ω) such that     Z Z Z   1   2 2 ¯ ¯ ¯    Ω µ (∇ × E)(∇ × F )dx − Ω (ω ² − jωσ)E · F dx − Ω (ω ² − jωσ)∇p · F dx = ¾ ½Z Z imp ¯  imp ¯  J · F dx + J · F dS ∀F ∈ HD (curl; Ω) −jω  S   Ω Γ2   ¾ ½Z Z Z     − (ω 2 ² − jωσ)E · ∇¯ J imp · ∇¯ q dx + J imp · ∇¯ q dS ∀q ∈ H 1 (Ω) . q dx = −jω  Ω

Ω

Γ2

S

D

(3.17) By repeating the trick with the substitution F = ∇q in the first equation, we learn that the Lagrange multiplier p identically vanishes, and for that reason, it is frequently called the hidden variable. In contrary to the original formulation, the stability constant for the regularized formulation converges to one, as k → 0. The regularized formulation works 1 because gradients of the scalar-valued potentials from HD (Ω) form precisely the null space of the curl-curl operator. The point about the regularized (mixed) formulation is that, whether we use it or not in the actual computations (the improved stability is one good reason to do it...), the original variational problem is equivalent to the mixed problem. This suggests that we cannot escape from the theory of mixed formulations when analyzing the problem. 11

3.3

Lagrange Multipliers

As we will see in this subsection, it turns out that the space of Lagrange multipliers (scalar potentials) becomes crucial for a number of interesting cases, such as steady state problems, axisymmetric problems, problems with harmonic dependence upon a spatial variable, etc. Variational formulation for steady state problems. If frequency ω = 0, then electric field E can be expressed as the gradient of a scalar potential, i.e., E = ∇p. ρ Thus, (2.4) becomes ∆p = . And the corresponding variational formulation is given by: ²  1 Find p ∈ HD (Ω) such that    Z Z 1 1   ρ¯ q dx ∀q ∈ HD (Ω) ∇p∇¯ q dx =  Ω

Ω

(3.18)

²

As we see, the problem can be easily solved in the space of scalar potentials. A similar situation occurs for steady state problems (2.5) and (2.6). Variational formulation for axisymmetric problems. Using system of equations (2.8), we can reduce the 3D variational formulation (3.11), to a 2D variational formulation for axisymmetric problems, given by:   Find (0, Eφ , 0) ∈ HD (curl; Ω) such that      Z   Z 1 ∂Eφ ∂ F¯φ 1 ∂(ρEφ ) ∂(ρF¯φ )        

Ω

µ

{

∂z ∂z

+

ρ2

= −jω

∂ρ

½Z

Ω

∂ρ

Jφimp F¯φ dx

}dx − +

Z

(ω ² − jωσ)Eφ F¯φ dx − jω 2

Ω

Γ2 ∪Γ3

imp ¯ JS,φ Fφ dS

¾

Z

Γ3

ˆ φ F¯φ dS β(ˆ n × φ)E

for all (0, Fφ , 0) ∈ HD (curl; Ω) . (3.19)

If domain Ω is convex, one can replace (0, Eφ , 0), (0, Fφ , 0) ∈ HD (curl; Ω) by Eφ , Fφ ∈ Thus, we can solve the problem in the space of scalar potentials.

1 HD (Ω).

Remark: If domain Ω is not convex then, in general, the above problem is not solvable in H 1 (Ω) (for details, see [1]). Nevertheless, one may solve for meridian components of magnetic field (Hρ , 0, Hz ) in H(curl(2D); Ω). Variational formulation for problems with harmonic dependence upon the spatial variable z. At this point, we assume that computational domain Ω is convex. Then, using system of equations (2.10), we can reduce the 3D variational formulation (3.11) to a 2.5D variational formulation given by: 12

 1  (Ω) such that Find (E 2D = (Ex , Ey ), p) ∈ HD (curl(2D); Ω) × HD     Z Z   1 1   ¯ 2D )dx + ¯ 2D dx  (∇ × E 2D ) · (∇ × F (jζE 2D − ∇p) · jζ F    Ω µ Ω µ   Z Z   2D  2D ¯ 2 ¯ 2D =   − (ω ² − jωσ)E · F dx − jω β(ˆ n × E 2D )F   Ω Γ3   ½Z ¾  Z   2D 2D

−jω

Ω

¯ J imp,2D · F

dx +

Γ2 ∪Γ3

¯ J imp,2D ·F S

dS

for all F ∈ HD (curl; Ω) .

   Z   1   (jζE 2D − ∇p) · (−∇¯ q )dx    µ Ω   Z Z    2   − (ω ² − jωσ)p¯ q dx − jω β(ˆ n × zˆ)¯ q dS =    Ω Γ3   ¾ ½Z Z    1   −jω JSimp,z · q¯dS for all q ∈ HD (Ω) . J imp,z · q¯dx + Ω

(3.20)

Γ2 ∪Γ3

1 (Ω). Again, we can take advantage of the space of scalar potentials, since p = Ez ∈ HD

13

4

A Fully Automatic hp-adaptive Finite Element (FE) Method

In order to solve a large variety of electromagnetic problems, we need a numerical technique that provides low discretization errors and, simultaneously, solves the discretized problem without prohibitive computational cost. In this context, an adaptive hp-Finite Element technique seems to satisfy both properties.

4.1

hp-Finite Elements (FE).

With each finite element, we associate element size h and order of approximation p. In the h-adaptive version of FE method, element size h may vary from element to element, while order of approximation p is fixed (usually p=1,2). In the p-adaptive version of FE method, p may vary locally, while h remains constant through the adaptive procedure. Finally, a true hp-adaptive version of FE method allows for both the varying of h and p locally. Two main reasons motivate the use of hp-FE for electromagnetic problems. 1. Lower order methods (e.g. h-FE) are suitable for rough (non-smooth) functions, but not for smooth functions. Conversely, higher order methods (e.g. p-FE) approximate very accurately smooth functions, but not rough functions. In real world applications, solutions usually exhibit both smooth and rough components. In this context, the flexibility of the hp-version of FE allows for approximating both parts of the solution very accurately. Indeed, the hp-version of FE guided with an adequate adaptive strategy shall provide exponential convergence rates for a large class of functions [28, 3, 2]. 2. The dispersion (pollution) error for wave propagation problems decreases considerably as p increases [18, 19]. Thus, the h-version of FE is not adequate for large frequencies. On the other hand, small elements are needed to capture geometrical details, such as antennas or thin edges. Only hp-FE combines both these features: small elements and high orders of approximation.

4.2

hp-edge Finite Elements. The de Rham diagram.

In order to solve Maxwell’s equations using finite elements, a family of H(curl) conforming elements is needed, since the space of admissible solutions is contained in H(curl). Nedelec introduced in 1980 [22] and 1986 [23] two families of FE that are conforming in the space H(curl;Ω)={E ∈ L2 (Ω) : ∇ × E ∈ L2 (Ω)}, Ω ⊂ R I 3 . These families were constructed 14

for both tetrahedrons and hexahedrons, and are referred to as Nedelec’s elements or edge elements (since for lowest order elements, degrees of freedom are associated to edges). In 1994, Monk [21] studied uniform hp-extensions of Nedelec’s curl conforming elements to establish interpolation error estimates for the hp-version of edge FE with constant p. A detailed description of edge elements for a true hp-FE code can be found in Demkowicz et al. [16, 6]. Crucial for the analysis of hp-edge elements is the so called de Rham diagram, which relates two exact sequences of spaces, on both continuous and discrete levels, and the corresponding interpolation operators. A short introduction to de Rham diagram for hp-edge elements can be found in the appendix of this report. hp-edge elements, resulting from the H(curl)-conforming FE approximation, involve the implementation of vector-valued shape functions with different dof, in order to take into account tangential and normal components of the electric field. Rachowicz and Demkowicz implemented both a 2D [26] and a 3D [27] hp-FE code for electromagnetics.

4.3

Survey of hp- and hp-edge FE codes.

A number of commercial and non-commercial hp-FE codes have been developed during the last two decades. • PHLEX, an advanced hp-adaptive finite element kernel trademark of COMCO Corp. It also includes an automatic adaptive procedure in both h and p, but decision between h- and p-refinement for a given element has to be performed manually. • PolyFEM, a trademark of Computer Aided Design Software, Corralville, Iowa. • Pro/MECHANICA, a trademark of Parametric Technology Corp., Waltham, Masachussets. • STRESSCHECK, a trademark of ESRD Res. Inc., St. Louis, Missouri. • STRIPE, developed by Flygtekniska F¨ors¨oksanstalden, Bromma, Sweden. • PZ, an hp-adaptive FE environment developed by Philippe Devloo, Unicamp, Brasil. • Philipp Frauenfelder and Ch. Schwab are developing an hp-FE code for general elliptic 3D problems, Z¨ urich, Switzerland. • 2Dhp90 and 3Dhp90 under development by Dr. Demkowicz and his team, Austin, Texas. 15

A 2D version of the last code will be used for this research. Next section lists some of the main features of this code.

4.4

2Dhp90 (3Dhp90): An anisotropic 2D (3D) fully automatic hp-adaptive package

The main features of the 2Dhp90 (3Dhp90) package are the following. • It utilizes isoparametric quads (hexahedras). • It allows for both isotropic and anisotropic mesh refinements. • It interfaces with a Geometrical Modeling Package (GMP) for an initial mesh generation and mesh refinements [31]. • It incorporates a new data structure in Fortran 90, described in [15]. • It reconstructs constrained element information by ascending and descending within trees of nodes, an abstraction for element vertexes, edges, and interiors. • It separates logical operations in the mesh refinement routines into two levels: 1. operations for nodes - problem independent. 2. operations for nodal degrees of freedom - problem dependent. • It incorporates a fully automatic hp-adaptive strategy [14]. The remainder of this section is devoted to a presentation on the integration of the two grid solver with the automatic hp-adaptive strategy. Indeed, the two grid solver is not intended to be used only as a solver itself, but also as a crucial part of the hp-adaptive strategy.

4.5

The fully automatic hp-adaptive strategy

As mentioned in the abstract, the hp-adaptive strategy presented in [14] produces a sequence of optimally hp-refined meshes that deliver exponential convergence rates in terms of size of the discrete problem (number of degrees-of-freedom (d.o.f.)). A given (coarse) hp-mesh is first refined globally both in h and p to yield a fine mesh, i.e. each element is broken into four element sons (eight in 3D), and the order of approximation is raised uniformly by one. Then, we solve the problem of interest on the fine mesh. The next optimal coarse mesh is then 16

determined by minimizing the projection based interpolation error of the fine mesh solution with respect to the optimally refined coarse mesh. The algorithm is very general, and it applies not only to H 1 -conforming but also to H(curl)- and H(div)-conforming discretizations as well [11, 13]. Moreover, since the mesh optimization process is based on minimizing the interpolation error rather than the residual, the algorithm is problem independent and can be applied to nonlinear and eigenvalue problems as well.

4.6

Goal-oriented hp-adaptivity

In the previous section, we outlined an algorithm that produces a sequence of optimal hpmeshes in the sense that each hp-grid minimizes the H 1 -norm of the error with respect the number of dof’s. But for a number of problems, as for example, modeling of some LWD electromagnetic measuring devices, we want to control concrete features of the solved problem (quantities of interest). These quantities can be represented, e.g., as bounded linear functionals of the solution. An adaptive algorithm based on a quantity of interest (different than the standard H 1 -norm) is called goal-oriented adaptivity, and it generates a sequence of meshes that is optimal in the sense that it minimizes the error of the quantity of interest. An automatic goal-oriented hp-adaptive algorithm has been implemented by Solin and Demkowicz [29], and it is similar to the hp-adaptive algorithm described in section 4.5. In particular, it also requires the solution of the fine grid problem.

4.7

The need for a two grid solver

Crucial to the success of both hp-adaptive strategies (the H 1 -norm adaptivity and the goaloriented adaptivity) is the solution of the fine grid problem. Typically in 3D, the global hp-refinement increases the problem size at least by one order of magnitude, making the use of an iterative solver inevitable. With a multigrid solver in mind, we implemented first a two-grid solver based on the interaction between the coarse and fine hp meshes. The choice is quite natural. The coarse meshes are minimum in size. Also, for wave propagation problems in the frequency domain, the size of the coarsest mesh in the multigrid algorithm is limited by the condition that the mesh has to resolve all eigenvalues below the frequency of interest. Consequently, the sequence of multigrid meshes may be limited to just a few meshes only.

17

5

Analysis of Edge Singularities Arising in Electrostatic Computations

A number of LWD tools incorporate electromagnetic (EM) antennas covered by metals with plastic apertures. Thus, edge singularities for the electric field may occur on the boundary between plastic and metal. As a result, to find the electric field near edge singularities may become essential. In addition, edge singularities for the electric (or magnetic) field may occur on the geological formation. Strength of an edge singularity is dependent upon geometry and sources. But the presented refinement algorithm do not only detects singularities, but also distinguish between different intensities of singularities. On this section, we focus on edge singularities arising in electrostatic problems, and we present high accuracy approximations for the electric field near this type of singularities. First, we consider a number of problems with edge singularities, for which analytical solution is known. These problems will be used as testing cases in order to asses performance of the adaptive strategy.

5.1

Model Problem

For our model problem, we consider two perfect electric conductor (PEC) infinite lines intersecting at a nonzero point (xc , yc ), as shown in Fig. 1. Let β be the angle between P EC1 and P EC2 , and ρ a Dirac’s delta function distribution of charges concentrated at the point (0, 0).

Boundary Value Problem (BVP). The electrostatic phenomena is governed by equation (2.4). Thus, solving for scalar potential p, we obtain: div∇p = ∆p =

ρ ²

With boundary conditions: • p = 0 on P EC1 and P EC2 . • p = 0 on a boundary far enough from the source. Since the electric field E decays as 1/r (where r is the distance from a given point to the source), for r large enough (for example, r = 106 ), the electric field magnitude is negligible, and can be replaced by 0. 18

C PE

2

(x c, yc )

PEC 1 Figure 1: Model problem Let Ω be the domain of interest (shown in Fig. 2). We denote by Γ to its boundary, i.e., ∂Ω = Γ. Then, our BVP is given by:

d=10 6

(−1,−1)

(−1,−1) 6

0 =1

d

d=10

6

β = 270 degrees

s

ee

gr

e 2d

d=2*10 6

β = 358 degrees

Figure 2: Computational domain Ω corresponding to angles β = 270 degrees and β = 358 degrees.   ∆p = ρ in Ω ²  p = 0 on Γ

19

(5.21)

Secondary field. Now, let ρ be a unitary electric charge distribution concentrated at (0,0) in free space. Solution is given by pprim = lnr. Thus, we obtain: (

∆p − pprim = 0 in Ω p − pprim = 0 − lnr = −lnr on Γ

(5.22)

Let psec = p − pprim . Then: (

∆psec = 0 in Ω psec = −lnr on Γ

(5.23)

By solving the problem for secondary potential psec , we avoid modelling of the source. Then, p = psec + pprim . In order to construct the variational formulation corresponding to the BVP (5.23), we multiply equation ∆psec = 0 by a test function v ∈ V = H01 (Ω) = {u ∈ H 1 (Ω) : u(∂Ω) = 0} and we integrate over the domain Ω. For the resulting integral to make sense, we introduce the space of admissible solutions U = {u ∈ H 1 (Ω) : u(∂Ω) = −lnr}. Finally, integration by parts yields to the following variational formulation problem:    u∈U Z   ∇u∇v = 0 Ω

∀v ∈ V

(5.24)

Exact Solution. The exact solution of this electrostatic problem can be computed analytically. At points located on the surface of P EC1 , the normal component of the electric field as a function of β, γ, x, y, xc , and yc , is given by (see [?]): π

l β sin( πγ ) 2π β En = − π s [1 − 2l β cos( πγ ) + l 2π β ] β

where s is the distance from a given point (x, y) to the corner (xc , yc ), i.e., s=

q

(x − xc )2 + (y − yc )2

and l is given by:

l=

q

x2c + yc2

20

s

In particular, for γ = β2 , we have: 1 2π −2π En = − q = π π π π π π −1− − −1+ 1+ β + l β] [(x2c + yc2 ) 2β s β + (x2c + yc2 ) 2β s1− β ] x2c + yc2 [l Notice that there exists a singularity at point (xc , yc ) (edge singularity in three dimensions) if and only if π < β. Furthermore, as larger is β as stronger is the singularity. In particular, strength of singularity is independent of the selected nonzero point (x c , yc ). Therefore, we set xc = yc = −1, obtaining: En = En (s, β) =

[2

−2π π π s + 2 2β s1− β ]

π 1+ π − 2β β

For simplicity, we will restrict ourselves to only two different angles β: 1. β = 270 degrees (physically interesting case). 2. β = 358 degrees (strong edge singularity).

5.2

Numerical Results Corresponding to Model Problem 1: 270 degrees

First, we present results corresponding to the model problem discussed above with β = 270 degrees. Convergence history can be found in Fig. 3. Straight line on the asymptotic regime reflects the exponential convergence of the method. Figures 4 and 5 show different zooms of the final hp-grid toward the singularity. Notice that elements close to the singularity are eleven orders of magnitude smaller than other elements. Finally, in figure 6, a comparison between the exact and approximate solution of the normal component of the electric field over P EC1 at points close to the singularity is presented. Since size of domain is 106 and we are focusing at points as close to the singularity as 10−6 , we need elements which are twelve order of magnitude smaller than the size of the domain. Since we stopped our algorithm when the smallest element was only eleven orders of magnitude smaller than the size of the domain, we do only obtain a good agreement between exact and approximate solutions up to a point 10−5 from the singularity.

21

2Dhp90: D A Fully automatic m hphp-adaptive Finite EElement m code 62.76error CALE

30.89 15.21 7.49 3.68 1.81 0.89 0.44 0.22 0.11 0.05 21

201

722

1763

3505

6126

9808

14730

nrdof 21074

Figure 3: Convergence history in the sca e nr. of dof (unknowns) (in the sca e dof to the power of 1/3) versus the H 1 -norm error in percentage (in the ogarithmic sca e)

5.3

Numerical Results Corresponding to Model Problem 2: 358 degrees

Now, we present results corresponding to the model problem discussed above with β = 358 degrees. Figures 7,8 and 9 show different zooms of the final hp-grid toward the singularity. Notice that elements close to the singularity are thirteen orders of magnitude smaller than other elements. Finally, in figure 10, a comparison between the exact and approximate solution of the normal component of the electric field over P EC1 at points close to the singularity is presented. Since size of domain is 106 and we are focusing at points as close to the singularity as 10−6 , we need elements which are twelve order of magnitude smaller than the size of the domain. Since we stopped our algorithm when the smallest element was thirteen orders of magnitude smaller than the size of the domain, we do fully resolve the problem on the region of interest.

22

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 1 D

A

m

Zoom = 10 m

D

A

Zoom = 102 D

A

m

m

m

Zoom = 103

m

D

Zoom = 104

A

m

m

Zoom = 105

Figure 4: Zooms of fina hp-grid toward the singu arity

23

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 106 D

A

m

Zoom = 107

m

D

A

Zoom = 108 D

A

m

m

m

Zoom = 109

m

D

Zoom = 1010

A

m

m

Zoom = 1011

Figure 5: Zooms of fina hp-grid toward the singu arity

24

40.491908

27.158549

39.656288

26.324481

38.820669

25.490413

37.985049

24.656345

37.149430

23.822277

36.313810

22.988209

35.478191

22.154141

34.642571

21.320073

33.806951

20.486005

32.971332

19.651937

32.135712

18.817869

31.300093

17.983801

30.464473

17.149733

29.628854

16.315665

28.793234

15.481597

27.957614 27.121995 -1.000001

14.647529 -1.000013

-1.000026

-1.000038

-1.000051

-1.000063

-1.000075

-1.000088

13.813461 -1.000100

x -1.000100

Value on the edge of normal component of electric field at distances 10−6 − 10−4 from the corner on the decibel scale.

-1.001338

-1.002575

-1.003813

-1.005050

-1.006288

-1.007525

-1.008763

x -1.010000

Value on the edge of normal component of electric field at distances 10−4 − 10−2 from the corner on the decibel scale.

13.831790 12.732826 11.633862 10.534898 9.435934 8.336970 7.238006 6.139042 5.040078 3.941114 2.842150 1.743186 0.644222 -0.454742 -1.553706 -2.652670 -3.751634 -1.010000

-1.133750

-1.257500

-1.381250

-1.505000

-1.628750

-1.752500

-1.876250

x -2.000000

Value on the edge of normal component of electric field at distances 10−2 − 100 from the corner on the decibel scale. Figure 6: Blue curve denotes exact solution, while FE approximation is represented by black curve.

6

Limitations of the fully automatic hp-adaptive strategy

On this section, we present some of the limitations of the fully automatic hp-adaptive strategy presented above, when we attempt to solve some real life Petroleum Engineering EM problems.

25

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 1 D

A

m

Zoom = 10 m

D

A

Zoom = 102 D

A

m

m

m

Zoom = 103

m

D

Zoom = 104

A

m

m

Zoom = 105

Figure 7: Zooms of fina hp-grid toward the singu arity

26

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 106 D

A

m

Zoom = 107

m

D

A

Zoom = 108 D

A

m

m

m

Zoom = 109

m

D

Zoom = 1010

A

m

m

Zoom = 1011

Figure 8: Zooms of fina hp-grid toward the singu arity

27

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 1012

Zoom = 1013

Figure 9: Zooms of fina hp-grid toward the singu arity

6.1

An axisymmetric battery in free space

Let H = (0, Hφ , 0) be the magnetic field produced by an axisymmetric battery in free space. And let Ω be the corresponding 2D domain, as shown in figure 11. Boundary Value Problem This EM problem is governed by Maxwell equations, which can be further reduced to the wave equation for the magnetic field: ∇ × σ −1 ∇ × H + jωµH = 0 where σ = σ + jω². Boundary conditions are given by:    Hφ = 1 on Γ1  

Hφ = 0 on Γ3 ∂Hφ = 0 on Γ2 ∪ Γ4 ∂n

(6.25)

Variational formulation Since domain Ω is not convex, Hφ may not belong to H 1 (see [1]). Therefore, H 1 -conforming FE shou d not be used as the space of admissib e so utions. However, one can prove that electric field E = (Eρ , 0, Ez ) ∈ H(curl(2D); Ω), and we can use edge FE (see [22, 23]) in order to find E. By doing so, our BVP in terms of the electric field becomes:

28

58.199630

38.328379

56.956577

37.080578

55.713524

35.832777

54.470471

34.584976

53.227418

33.337175

51.984366

32.089374

50.741313

30.841572

49.498260

29.593771

48.255207

28.345970

47.012154

27.098169

45.769102

25.850368

44.526049

24.602567

43.282996

23.354765

42.039943

22.106964

40.796891

20.859163

39.553838 38.310785 -1.000001

19.611362 -1.000010

-1.000019

-1.000027

-1.000036

-1.000045

-1.000054

-1.000063

18.363561 -1.000072

x -1.000072

Value on the edge of normal component of electric field at distances 10−6 − 10−4 from the corner on the decibel scale.

-1.000962

-1.001852

-1.002742

-1.003633

-1.004523

-1.005413

-1.006303

x -1.007193

Value on the edge of normal component of electric field at distances 10−4 − 10−2 from the corner on the decibel scale.

18.385059 16.854529 15.323999 13.793469 12.262939 10.732409 9.201879 7.671349 6.140819 4.610289 3.079759 1.549229 0.018699 -1.511831 -3.042361 -4.572891 -6.103421 -1.007193

-1.096212

-1.185230

-1.274248

-1.363267

-1.452285

-1.541303

-1.630322

x -1.719340

Value on the edge of normal component of electric field at distances 10−2 − 100 from the corner on the decibel scale. Figure 10: Blue curve denotes exact solution, while FE approximation is represented by black curve.   ∇ × ( µ1 ∇ × E) − (ω 2 ² − jωσ)E = 0    

n × ∇ × E = −jωµ on Γ1

 n × ∇ × E = 0 on Γ3     n × E = 0 on Γ ∪ Γ 2 4

And the corresponding variational formulation is given by:

29

(6.26)

4

3

1

d=4 m

d=0.01 m

2

2 d=0.01 m

4

3 z

d=2 m

4 r

Figure 11: A domain Ω representing a battery in free space

  Find E ∈ HD (curl(2D); Ω) such that      Z Z   1 2        

Ω

µ

¯ )dx − (∇ × E) · (∇ × F ω2µ

½Z

¯ dS F

Γ1

Ω

¯ dx = (ω ² − jωσ)E · F

¾

for all F ∈ HD (curl(2D); Ω) .

where HD (curl(2D); Ω) = {E ∈ H(curl; Ω) : n × E(Γ2 ∪ Γ4 ) = 0} is the space of admissible solutions.

6.2

hp-edge Finite Element approximations using the fully automatic hp-adaptive strategy

At this point, we discretize the problem using hp-edge FE, as described in [6, 27]. First, we construct a rather coarse initial grid (see figure 12). Then, we execute the fully automatic hp-adaptive algorithm described on previous sections. As mentioned before, the objective of this algorithm is to minimize the energy norm of the error (with respect the number of unknowns). And as expected, exponential convergence rates are obtained (Fig. 13). 30

2Dhp90: A Fu 2Dhp90 Fullyy automatic u om hp hp-adaptive d p v FFinite n E Element m n code od

y z x

Figure 12: A zoom of the initia hp-grid.

2Dhp90 A Fu y au oma c hp adap ve F n e E emen code 90 19e o SCALES n do ^0 33 og e o

71 92 57 35 45 74 36 47 29 09 23 20 18 50 14 75 11 76 9 38 27

183

579

1325

2534

4317

6786

Figure 13: Convergence history.

31

10053

n do 14229

By looking at the final hp-grid (fig. 14,15), we observe a large number of refinements taking place around the two singular points of the domain, which are located at the vertexes of the battery. Indeed, in order to minimize the energy norm of the error, it is required to approximate singularities rather accurately. Unfortunately, the main objective of this type of problems is to estimate the value of the magnetic field at points far from the battery. Thus, the final hp-grid constructed using the fully automatic hp-adaptive algorithm is inadequate for our pourposes, and it becomes useless. A fully automatic goal-oriented hp-adaptive algorithm is needed.

32

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 101

Zoom = 1 D

A

m

m

D

A

Zoom = 102 D

A

m

m

m

Zoom = 103

m

D

Zoom = 104

A

m

m

Zoom = 105

Figure 14: Zooms of fina hp-grid toward a singu ar point.

33

2Dhp90: D A Fully automatic m hp-adaptive Fin Finite Element m code

D

A

m

m

y z x

Zoom = 106 D

A

m

Zoom = 107

m

D

Zoom = 108

A

m

m

Zoom = 109

Figure 15: Zooms of fina hp-grid towards a singu ar point.

7

Summary and Conclusions

On this report, we first discussed Maxwell’s equations in the frequency domain, with a particular emphasis on a number of Petroleum Engineering EM applications. Then, the corresponding variational formulation was derived. It turns out that the space of scalar potentials (Lagrange multipliers) becomes crucial for a variety of interesting EM problems. Thus, H 1 and/or H(curl) spaces are required for solving EM problems. hp- discretizations of H 1 and H(curl) spaces were discussed in Section 4. At that point, we presented the fully automatic hp-adaptive strategy, which minimizes the energy norm of the error by automatically generating a sequence of optimally refined grids. This sequence of grids provides exponential convergence rates with respect to the number of unknowns (d.o.f). 34

Then, we utilized the adaptive strategy in order to solve an electrostatic model problem with edge singularities of different strength. We proved that high accuracy solution of such a problem is possible due to the use of higher order finite element methods combined with a fully automatic adaptive procedure. This numerical technique has already been successfully implemented for electromagnetic problems as well as for three dimensional problems. Unfortunately, as observed in Section 6, for a number of Petroleum Engineering EM applications, the adaptive procedure should not be oriented toward minimization of the energy norm of the error, but toward minimization of an arbitrary quantity of interest (goal-oriented adaptivity). In this sense, it is desirable to implement a fully automatic hp-goal oriented adaptive algorithm based on the energy driven algorithm already implemented.

8

Future Work

In this section, I briefly discuss the required work that, in my opinion, is needed in order to transform the current hp-adaptive FE code into a powerful tool for solving Petroleum Engineering EM problems. I am firmly convinced that a large number of real life Petroleum Engineering EM problems can be solved with an optimal grid (in goal) produced by an hp goal-oriented adaptive algorithm. From the mathematical point of view, goal-oriented adaptivity is already well understood. From the software development point of view, a first version for H 1 problems was already implemented by Dr. Demkowicz et al. in 2002 [29]. But since a new and improved version of the energy driven adaptive algorithm was released on December 2002, it would be desirable to produce also a new version of the goal-oriented adaptive algorithm. This new version should incorporate discretization of both H 1 and H(curl) spaces. As shown in [29], the algorithm shall certainly produce a sequence of optimal grids with respect an arbitrary goal (target) introduced by the user. But crucial to the success of both energy driven and goal-oriented hp-adaptive strategies is the ability of solving the fine grid problem. As mentioned on section 4, in 3D typically, the number of unknowns of the fine grid problem is at least one order of magnitude larger than the number of unknowns of the corresponding coarse grid problem. Also, notice that the average number of nonzero entries per column of an hp-grid is about twenty times larger than the number of nonzero entries per column of the corresponding h-grid. In this context, a two grid solver becomes a necessity rather than an option. A 2D and 3D hp two grid solver for H 1 problems has already been implemented (see 35

[25] for details), while there is an ongoing effort to produce a similar algorithm for H(curl problems. Finally, the current implementation of the geometrical modeling package (GMP) [31] already allows for description of rather complex 2D and 3D logging tools. In my opinion, no improvement is necessary in that area. In terms of graphics, visualization of a large variety of 1D and 2D graphics as well as a number of 3D graphics is currently available. Final Thought: In my opinion, a qualitative improvement in solving complex Petroleum Engineering EM problems will come from combining the expertise froma number of different areas of sciences, such as Mathematics, Numerical Analysis, Geophysics, Engineering, etc.

36

9

Appendix: de Rham diagram

In this appendix, we define formally a number of (scalar-valued and vector-valued) polynomial spaces used for implementation of a variational formulation in terms of hp- and hp-edge finite elements. Then, we define three interpolation operators from continuous spaces (defined below) to corresponding discretized polynomial subspaces. Finally, we present de Rham diagram. This appendix is intended only to introduce de Rham diagram, and details are not discussed. For a more rigorous and detailed presentation, see [11, 12]. Let K = (0, 1)3 be the master cube. We introduce the following polynomial spaces over K. • Qp = Q(p

1 ,p2 ,p3 )

- space of polynomials over K of order pi with respect to xi .

• W p = Qppf ,pe - space of polynomials of Qp , whose traces to faces f reduce to (possibly lower) order pf = (p1f , p2f ), f = 1, . . . , 6, and traces to edges e reduce to (possibly lower) order pe , e = 1, . . . , 12. Note that (possible lower) order is meant with respect the polynomial degree on the parallel direction to the face/edge. This condition is implicitly understood in all subspaces defined below. p p = Qp−1 • W0,K = Qp−1,−1 ; W0,f

• Qp = {E ∈ Qp+(0,1,1) ×Qp+(1,0,1) ×Qp+(1,1,0) , Et |e ∈ P pe (e), Et |f ∈ Qpf +(0,1) ×Qpf +(1,0) } • Qp0,K = {E ∈ Qp+(0,1,1) × Qp+(1,0,1) × Qp+(1,1,0) , Et |e = 0, Et |f = 0} • V p = {u ∈ Qp+(1,0,0) × Qp+(0,1,0) × Qp+(0,0,1) , un |f ∈ Qpf } Now, we define three interpolation operators. Π : W = H 1/2+r (K)(r > 0) −→ Qp where Πu = uhp is defined as follows.  uhp (v) = u(v)          |uhp − u|0,e → min

   |uhp − u| 1 ,f → min   2    

|uhp − u|1,K → min 1

for each vertex v for each edge e for each face f for cube K ,

Πcurl : Q = {E : E ∈ H 2 +r (K), ∇ × E ∈ H r (K), r > 0} −→ Qp 37

where Πcurl E = E hp is defined as follows.  |Ehp − E|−1,e → min         |(∇ × E hp ) · nf − (∇ × E) · nf |− 1 ,f → min   2   

(E

− E, ∇ φ)

1

∀φ ∈ W p

hp f − 2 ,f 0,f        |∇ × E hp − ∇ × E|0,K → min       p

(E hp − E, ∇φ)0,K

∀φ ∈ W0,K

for each edge e

for each face f

for cube K .

Πdiv : V = {u : u ∈ L2 (K), ∇ ◦ u ∈ H r (K), r > 0} −→ V p where Πdiv u = uhp is defined as follows.  kuhp · nf − u · nf k− 1 ,f → min   2   

for each face f

k∇ ◦ uhp − ∇ ◦ uk0,K → min

    

∀φ ∈ Qp0,K

(uhp − u, ∇ × φ)0,K

for cube K .

Finally, we present de Rham diagram (in 3D): R I −→   yid

W   yΠ

∇ −→

Q

∇× −→

V

∇◦ −→

  div yΠ

  curl yΠ

L2   yP

−→

0

∇× ∇ ∇◦ R I −→ W p −→ Qp −→ Vp −→ W p−1 −→ 0 .

38

(9.27)

References [1] F. Assous, P. Ciarlet, Jr., S. Labrunie, ”Theoretical tools to solve the axisymmetric Maxwell equations”, Math. Meth. Appl. Sci., Vol 25, pp. 49-78, (2002). [2] I. Babuˇska, B. Guo, ”Approximation properties of the hp-version of the finite element method”, Comp. Meth. Appl. Mech. Engg., 133 (1996), 319-346. [3] I. Babuˇska, B. Guo, ”Regularity of the Solutions of Elliptic Problems with Piecewise Analytic Data, parts I and II”, SIAM J. Math. Anal. 19 (1988), 172-203 and 20, (1989), 763-781. [4] C. Balanis, ”Advanced Engineering Electromagnetics”, John Wiley & Sons, 1989. [5] M. Cessenat, Mathematical Methods in Electromagnetism, Linear Theory and Applications, World Scientific, London 1996. [6] L. Demkowicz, ”Edge Finite Elements of Variable Order for Maxwell’s Equations: A Discussion”, TICAM Report 00-32, 2000. [7] L. Demkowicz, ”2D hp-Adaptive Finite Element Package (2Dhp90). Version 2.0”, TICAM Report 02-06, January 02. [8] L. Demkowicz, ”Finite Element Methods for Maxwell Equations”, In preparation for publication in Encyclopedia of Computational Mechanics, 2003. [9] L. Demkowicz, J.T. Oden, W. Rachowicz, O. Hardy, ”Toward a Universal h−p Adaptive Finite Element Strategy, Part 1. Constrained Approximation and Data Structure.”, Texas Institute for Computational Mechanics, The University of Texas at Austin, March 1989. [10] L. Demkowicz, K. Gerdes, Ch. Schwab, A. Bajer, and T. Walsh, ”HP90: A General and Flexible Fortran 90 hp-FE Code”, Computing and Visualization in Science, 1, 145-163, 1998. [11] L. Demkowicz, P. Monk, L. Vardapetyan, W. Rachowicz, ”De Rham Diagram for hp Finite Element Spaces”, Ticam Report 99-06, 1999 [12] L. Demkowicz, I. Babuˇska, J. Schoberl, P. Monk, ”De Rham Diagram in 3D. Quasi Optimal p-Interpolation Estimates”, Ticam Report, 2002 [13] L. Demkowicz, and I. Babuˇska, ”Optimal p Interpolation Error Estimates for Edge Finite Elements of Variable Order in 2D”, TICAM Report 01-11, submitted to SIAM Journal on Numerical Analysis 39

[14] L. Demkowicz, W. Rachowicz, P. Devloo, ”A Fully Automatic hp-Adaptivity”, Journal of Scientific Computing, 17, 1-3, 127-155, 2002. [15] L. Demkowicz, D. Pardo, W. Rachowicz, ”3D hp-Adaptive Finite Element Package (3Dhp90). Version 2.0. The Ultimate (?) Data Structure for Three-Dimensional, Anisotropic hp-Refinements”, Ticam Report 02-24, June 2002 [16] L. Demkowicz, L. Vardapetyan, ”Modeling of Electromagnetic Absortion/Scattering Problems Using hp-adaptive Finite Elements”, Comput. Methods Appl. Mech. Engrg., 152 (1998) 103-124. [17] W. Gui, I. Babuˇska, ”The h, p and h − p Versions of the Finite Element Method in One Dimension. Part 1. The Error Analysis of the p-Version. Part 2. The Error Analysis of the h- and h − p Versions. Part 3. The Adaptive h − p Version.”, Numer. Math. 49 (1986), 577-612, 613-657, 659-683. [18] F. Ihlenburg, I. Babuˇska, ”Finite Element Solution to the Helmoltz Equation with High Wave Numbers - Part I: the h-version of the FEM”, Computers Math. Applic. Vol 30, No. 9, 9-37, 1995. [19] F. Ihlenburg, I. Babuˇska, ”Finite Element Solution to the Helmoltz Equation with High Wave Numbers - Part II: the h − p-version of the FEM”, SIAM J. Numer. Anal. Vol 34, No. 1, 315-358, 1997. [20] J. R. Lovell ”Finite Element Methods in Resistivity Logging”, Deft University of Technology, 38-45,1993. [21] P. Monk, ”On the p- and hp-extension of Nedelec’s curl-conforming elements”, J. Comput. Appl. Math. 53 (1994) 117-137. [22] J.C. Nedelec, ”Mixed Finite Elements in R I 3 ”, Numer. Math., 35, 315-341, 1980. [23] J.C. Nedelec, “A New Family of Mixed Finite Elements in R I 3 ”, Numer. Math., 50, 57-81, 1986. [24] J.T. Oden, L. Demkowicz, W. Rachowicz, T.A. Westermann, ”Toward a Universal h − p Adaptive Finite Element Strategy, Part 2. A Posteriori Error.”, Texas Institute for Computational Mechanics, The University of Texas at Austin, March 1989. [25] D. Pardo, L. Demkowicz ”Integration of hp-Adaptivity and Multigrid. I. A Two Grid Solver for hp Finite Elements”, Ticam Report 02-33, Texas, 2002.

40

[26] W. Rachowicz, L. Demkowicz, ”A Two-Dimensional hp-Adaptive Finite Element Package for Electromagnetics (2Dhp90 EM)”, Ticam Report 98-16. [27] W. Rachowicz, L. Demkowicz, ”An hp-Adaptive Finite Element Method for Electromagnetics”, Int. J. Numer. Meth. Engng 2002; 53: 147-180. [28] Ch. Schwab, ”p- and hp- Finite Element Methods”, Oxford Science Publications. [29] P. Solin, L. Demkowicz Goal-Oriented hp-Adaptivity for Elliptic Problems, Ticam Report 02-32, August 2002. [30] A. Taflove, ”Why Study Electromagnetics: The First Unit in an Undergraduate Electromagnetics Course”, Department of Electrical and Computer Engineering, Northwestern University, Evanston, Il, 2002. [31] D. Xue, L. Demkowicz ”Geometrical Modeling Package. Version 2.0”, TICAM Report 02-30,2002. [32] A. Zdunek, W. Rachowicz, ”A Three-dimensional hp-Adaptive Finite Element Approach to Radar Scattering Problems”, Presented ad Fifth World Congress on Computational Mechanics, Vienna, Austria, 7-12 July, 2002.

41