1 Introduction
New computational techniques meet the needs of industry in developing efficient computational methods to simulate partial differential equations (PDEs). Especially, for simulations in higher dimensions and large computational domains. In this regard, certain
domain decomposition (DD) method leading to efficient schemes, in numerical investigations, has gained a lot of interest in numerical analysis community. This variant of the DD method is the subject of our current, and some related and ongoing, research. This type of schemes was previously studied, e.g. in [12, 30] for others than studied here problems.The present work is a further development of the DD hybrid finite element/finite difference (FE/FD) method for timedependent Maxwell’s equations for electric field in nonconductive media, studied in [3, 4]. A stable, timedomain (TD), DDM scheme for Maxwell’s equations was proposed in [26, 27], and further verified in [13]. This method uses FDTD scheme on the, structured, FD part of the mesh, and edge elements on unstructured part. In applications, because of edge elements implementation, this method remains computationally expensive.
The domain decomposition FE/FD method for timedependent Maxwell’s equations for electric field, assuming constant dielectric permittivity function in a finite difference domain, was considered in [3]
. This assumption simplifies the numerical schemes in both FE and FD domains and significantly reduces computational efforts for implementation of the whole DD method. Modified numerical scheme, energy estimate and numerical verifications of this method was presented in
[4]. However, the fully stability and convergence analysis with numerical implementations, in  and norms of the developed FE and FD schemes, are not presented in the above studies. We fill this gap in the present work.More specificaly, we present stability analysis for explicit schemes for both FEM and FDM in the DD hybrid FE/FD method. The DDM is constructed such that FEM and FDM coincide on the common, structured, overlapping layer between the two subdomains. The resulting domain decomposition approach at the overlapping layers can be viewed as a FE scheme which avoids instabilities at the interfaces. Similar to the DD approach of [4, 5], we decompose the computational domain such that FEM and FDM are used in different subdomains: FDM in simple geometry and FE in the subdomain where more detailed information is needed about the structure of this subdomain. This also allows application of adaptive FEM in such subdomain, see, e.g. [2, 3, 8, 9, 10, 28, 29].
Reliability and convergence of the domain decomposition method, studied in this work, are evident for solution of coefficient inverse problems (CIPs) in , see, e.g. [8, 9, 10, 28, 29]. For the case of CIPs, the computational domain is splitted into subdomains such that a simple discretization scheme can be used in a large region and more refined discretization scheme is applied in smaller, however more critical, part of the domain. In most algorithms for solution of electromagnetic CIPs, to determine the dielectric permittivity function inside a computational domain, a qualitative collection of experimental measurements is necessary on it’s boundary or in it’s neighborhood. In such cases it is convenient to condsider the numerical solution of timedependent Maxwell’s equations in different subdomains with constant dielectric permittivity function in some subdomain and nonconstant in the other ones. For the timedependent Maxwell’s equations, the DD scheme of [4], which is analyzed in the present work, is used for solution of different CIPs to determine the dielectric permittivity function in nonconductive media using simulated and experimentally generated data, see [2, 8, 9, 10, 28, 29].
An outline of this paper is as follows. In Section 2 we introduce the mathematical model. In Section 3 we briefly present the domain decomposition FE/FD method and communication scheme between two methods. In Section 4 we describe the domain decomposition FE/FD method for solution of Maxwell’s equations and set up the finite element and finite difference schemes. Section 5 is devoted to the stability analysis. In Section 6, we derive optimal a priori error estimates in finite element method for the semidiscrete (spatial discretizations) problems. Finally, in Section 7 we present numerical implementations that justify the theoretical investigations of the paper. In what follows, will be a generic constant independent of all parameters, unless otherwise specifically specified, and not necessarliy the same at each occurance.
2 The mathematical model
The Cauchy problem for the electric field , , , of the Maxwell’s equations, under the assumptions that the dimensionless relative magnetic permeability of the medium is and the electric volume charges are zero, is given by
(1) 
where, and are the dimensionless relative dielectric permittivity and electric conductivity functions, respectively. , and are the permittivity and permeability of the free space, respectively, and is the speed of light in free space. In this paper we consider the problem (1) in nonconductive media, i.e. , and hence study the initial value problem
(2) 
To solve the problem (2) numerically, we consider it in a bounded domain (instead of whole ), with boundary , and employ a split scheme on : a hybrid, finite element/finite difference scheme, kind of domain decomposition, developed in [3, 4] and summarized in Algorithm 1. More specifically, we divide the computational domain into two subregions, and such that and is a subset of the convex hul of . The function is assumed to be constant in , and bounded and smooth in . The communication ß between and is arranged using an overlapping mesh structure through a twoelement thick layer around as shown by blue and green common boundaries in Figure 1. The blue boundary is outer boundary of and inner boundary of . Similarly, the green boundary is the inner boundary of from which the solution is copied to the green boundary of .
The key idea with such a decomposition is to be able apply different numerical methods in different computational domains. For the numerical solution of (2) in we use the finite difference method on a structured mesh. In , we use finite elements on a sequence of unstructured meshes , with elements consisting of triangles in and tetrahedra in , both satisfying minimal angle condition. This approach combines the flexibility of the finite elements and the efficiency of the finite differences in terms of speed and memory usage and fits well for reconstruction algorithms presented below.
We assume that for some known constant , the function satisfies
(3) 
Conditions (3) on and the relation
(4) 
together with divergence free field , make equations in (2) independent of each others in and so that, in , we just need to solve the system of wave equations:
(5) 
Remark
It is well known that, for stable implementation of the finite element solution of Maxwell’s equations, divergencefree edge elements are the most satisfactory ones from a theoretical point of view [21, 24]. However, the edge elements are less attractive for solution of timedependent problems, since a linear system of equations should be solved at each time iteration step. In contrary, P1 elements can be efficiently used in a fully explicit finite element scheme with lumped mass matrix [11, 14, 18]. It is also well known that numerical solution of Maxwell’s equations using nodal finite elements is often unstable and results spurious oscillatory solutions [22, 25]. There are a number of techniques to overcome such instabilities, see, e.g. [15, 16, 17, 23, 25].
In [6, 7], a finite element analysis shows stability and consistency of the stabilized finite element method for the solution of (1) with . In the current study we show stability and convergence for the combined FEM/FDM scheme, under the condition (3) on , where the stabilized FEM is used for the numerical solution of (2) in and usual FDM discretization of (5) is applied in .
Remark
Here, we consider the case when for . Further, we assume that , and . Recall that we assumed nonconductive media: . In the presence of electric conductivity, additional terms appear in the equations. They lead to more involved estimates and heavier implementations which we plan to perform in a forthcoming study.
Hence, in this note we study the following initial boundary value problem:
(6) 
3 The structure of domain decomposition
a)  b)  c) 
We now describe the DD method between two domains, and , where the FEM is used for computation of the solution in , and FDM is used in . Communication between and is achieved letting overlapping of both meshes across a twoelement thick layer around  see Figure 1. The common nodes of both and domains belong to either of the following boundaries (see Figure 1):

Nodes on the blue boundary  lie on the boundary of and are interior to ,

Nodes on the green boundary  lie on the inner boundary of and are interior to .
Then the main loop in time for the explicit hybrid FEM/FDM scheme, that solves (2) associates with appropriate boundary conditions, at each time step is described in Algorithm 1 below:
By (3), at the overlapping nodes between and . Thus, FEM and FDM schemes coincide on the common, structured, overlapping layer. Hence, we avoid instabilities at interfaces.
4 Derivation of computational schemes
In this section we construct, combined, finite elementfinite difference schemes to solve the model problem (6). To do this, first we present the finite element scheme to solve (6) in entire . This induces the finite element scheme in . Then, we derive the finite difference scheme in , when the domain decomposition FE/FD structure is applied to solve (6) in .
Remark
The computational schemes derived in this section are explicit and therefore, for their converegence, the CFL condition below (see ,e.g. [6, 7]) should be satisfied
(7) 
Here is a mesh independent constant, is the time step, and is the mesh size.
In the sequel, we denote the inner product of by , and the corresponding norm by . The scalar inner product in we denote by , and the associated norm by . Further, we let to be the boundary of , the inner boundary of and the outer boundary of .
4.1 Finite element discretization in
First we derive finite element scheme to solve the model problem (6) in whole . Next, we discretize in two steps: (i) the spatial discretization using a partition of into elements , where is a mesh function defined as , representing the local diameter of elements. We also denote by a partition of the boundary into boundaries of the elements such that, at least one of the vertices of these elements belong to . (ii) As for temporal discretization, we let be a uniform partition of the time interval into subintervals of length As usual, we also assume a minimal angle condition on elements in . To formulate the finite element method for (6) in , we introduce the finite element space for each component of the electric field defined by
where denote the set of piecewiselinear functions on . Setting we define (resp. ) to be the usual
interpolant of
(resp. ) in (6). Also, because of the Dirichlet boundary data in (6) we need to choose the test function space as(8) 
Then, recalling (4), and the fact that the spatial semidiscrete problem in reads:
Find such that ,
(9) 
We note that
implies that
(10) 
Recalling (8), vanishing test functions at the boundary yields , and the final weak formulation for the semidiscrete problem in is: Find such that ,
(11) 
For the, reflexive, inhomogeneous boundary condition, see the FE scheme for .
To get fully discrete scheme for (6) we apply time discretization to (9) approximating , denoted by , where we use the central difference scheme, for :
(12) 
Multiplying both sides of (12) by we get, that
(13) 
Rearranging terms in (13) we get the following scheme: Given the approximate initial data, and , find such that ,
(14) 
4.2 Finite element discretization in
To solve the model problem (6) via the domain decomposition FE/FD method, we use the split , see Figure 1. Thus in we use FEM to solve the equation
(15) 
Here, is the restriction of the solution obtained by the FDM in to and therefore the test functions are not vanishing at the boundary and hence the term corresponding to in (9) will be appearing in the weak formulation.
To formulate the finite element method for (15) in , mimiking (9), we introduce the finite element space for each component of the electric field defined by
Setting we define , , and to be the usual interpolants of , , and , respectively, in . Then, similar to the FE scheme for , we get the following finite element scheme for : Given , , and , find such that
(16) 
A corresponding fully discrete problem in reads as follows:
Given , , , , and ; find such that
(17) 
Remark
Note that, in (15), Dirichlet boundary condition can be considered as well.
4.3 Fully discrete FE scheme for the electric field in
We expand the functions in terms of the standard continuous piecewise linear functions in space as
(18) 
where denote unknown coefficients at and the spatial mesh point . Then, substituting of (18) in (14), and setting , we obtain the linear system of equations:
(19) 
Note that, unlike (14), now the contributions at boundary of the element appear in . Here, are the block mass matrices in space, are the block stiffness matrices in space, denote the nodal values of , and is a uniform time step. Now we define the mapping from the reference element onto such that and let be the piecewise linear local basis function on such that . Then, the explicit formulas for the entries in system of equations (19), for each element , can be written as:
(20) 
where , and , denote the , and , scalar products on and , respectively. Note that here is only the part of the boundary of element that lies at .
4.4 Fully discrete scheme for the electric field in
As in the fully discrete FE scheme (19) in , we obtain fully discrete FE scheme in in the domain decomposition setting: Expanding the functions of via the continuous piecewise linear functions in space as in (18), and then substituting them in (17), (with , and ), we get the linear system of equations:
(22) 
Here, is the block mass matrices in space, restricted to , otherwise the same as in (19), are the block stiffness matrices in space as in (19),
is the assembled load vector ,
denote the nodal values of , is the time step. All quantities are for . Defining the mapping for the reference element in the mesh generated in as in the previous section, the formulas for entries of all matrices in the system (22) are the same as those in (20), and the entries of load vector are computed as(23) 
Again, approximating with the lumped mass matrix , we obtain the following fully explicit scheme:
(24) 
4.5 Finite difference formulation
We recall now that from conditions (3) it follows that in the function . This means that in for the model problem (2) the forward problem will be
(25)  
(26)  
(27)  
(28) 
Using standard finite difference discretization of the equation (25) in we obtain the following explicit scheme for the solution of the forward problem:
(29) 
In the, system of, equations above, is the solution at the time iteration at the discrete point , is the time step, and is the discrete Laplacian.
Note that, in (29), the Dirichlet boundary consitions can be considered as well.
5 Stability
In this section we derive stability estimates for the semidiscrete approximations. For stability in these estimates are extensions of the stability approach derived for the continuous problem in [4]. As for the stability in we get slightly different norms involving contributions corresponding to the reflexive boundary: . We use discrete version of a triple norm induced by the weak variational formulation of (6), where we use the relation (10) (which is not necessary in the continuous case where , however, in general ):
Find such that
(30) 
Remark
In general, in nondivergent free case, the bilinear form induced by (30) is not coercive. Further conforming finite element may result in spurious solutions. A remedy is through modifying the equation by adding a gauge constrain of Coulombtype, see, e.g. [21] and [23]. This is supplied by the ”zero”term: , in (6), which we add in the continuous variational formulation in (9). This, however, is not necessarily true in the discrete forms, e.g. in (16), where most likely Taking in (30), (we used the boundary condition on ), yields
(31) 
which, due to the fact that is independent of , can be rewritten as
(32) 
Proposition [4]
Let be a bounded domain with piecewise linear boundary . Then, the equation (6) has a unique solution . Further Let and , then there is a constant such that, , the following stability estimate holds true
(33) 
Proof.
Below we translate this stability to the semidiscrete problem.
5.1 Stability estimate for the semidiscrete problem in
The stability for the semidiscrete problem in is basically as in the continuous case above where all :s are replaced by with some relevant assumptions in the discrete data, viz.
Lemma
Assume that the interpolants of the data and : and satisfy the regularity conditions and , then for each
(34) 
where
(35) 
5.2 Stability of the semidiscrete problem in
The stability of the semidiscrete problem in , relying on the variational formulation (16), and due to the appearance of the data function , is slightly different from (34). We rewrite (16), in view of (10), and with as: Given , , and , find such that
(36) 
To deal with the term we rewrite (36) in its original form as (9) for and with :
Comments
There are no comments yet.