Local FE approximation in Dirichlet DD methods

by Yi Zhang

After decomposing a domain, we can formulate approximations by FE or BE methods in local subdomains. In the case of Dirichlet domain decomposition method, the global unknown on the skeleton is primal, i.e. u(x), while local unknowns are the nodal variable within subdomain \Omega_i (FE) or on the local boundaries \Gamma_i(BE). If focusing on a specific subdomain, one can neglect the subscripts and assume model problem as

-\nabla^2u(x)=f,\quad x\in\Omega

u(x)=g_D,\quad x\in\Gamma_D

\partial u/\partial n=g_N,\quad x\in\Gamma_N

By using Green’s identity, the weak formulation of this problem is

a(u,v)-\int_{\Gamma}{v\nabla u\cdot}ds=\int_{\Omega}{fv}d\Omega=\langle f,v\rangle

In which a(u,v) is the bilinear form for space H^1_0(\Omega). Assigning test function satisfying v|_{\Gamma_D}=0, we have

a(u,v)=\langle f,v \rangle+\langle g_N,v \rangle_{\Gamma_N}

Now we decompose the primal u(x) into three components

u=u_0+P\tilde{u}_1+P\tilde{g}_D, \quad u_0\in H^1_0(\Omega),\tilde{u}_1\in H^{1/2}_0(\Gamma,\Gamma_D),\tilde{g}_D|_{\Gamma_D}=g_D

which means, the unknown variable consists of its value within the domain, extension/prolongation from unknown Dirichlet boundary value on \Gamma_N, and extension/prolongation from Dirichlet data.  Here P denotes the “extension/prolongation” operator, mapping a function as boundary value to a function on the domain. By this decomposition one has two unknowns, u_0 on the domain and \tilde{u}_1 on the boundary, so another weak formulation is needed. Assuming u_0 is the weak solution of governing equation, it is easy to show, for test function space H^1_0(\Omega)

a(u_0+P\tilde{u}_1+P\tilde{g}_D,v)=\langle f,v \rangle

Finally, the weak formulation of the model problem is to resolve couple equations

a(u_0+P\tilde{u}_1,v)=\langle f,v \rangle-\langle P\tilde{g}_D,v \rangle,\quad \forall v \in H^1_0(\Omega)

a(u_0+P\tilde{u}_1,Pw)=\langle f,Pw\rangle+\langle g_N,w \rangle_{\Gamma_N}-\langle P\tilde{g}_D,Pw\rangle,\quad \forall w\in H^{1/2}_0(\Gamma,\Gamma_D)

for unknown u_0 and \tilde{u}_1. Assigning two finite element spaces, one on \Omega, one on \Gamma:

X=\text{span}(\phi_i)\subset H^1_0(\Omega),\quad \tilde{X}=\text{span}(\psi_j)\subset H^{1/2}_0(\Gamma,\Gamma_D)

The Galerkin FE formulation of the coupled problem is then

a(\phi_i,\phi_j)u^i_0+a(P\psi_i,\phi_j)\tilde{u}_1^i=\langle f,\phi_j \rangle-\langle P\tilde{g}_D,\phi_j \rangle, \forall \phi_j\in X

a(\phi_i,P\psi_j)u_0^i+a(P\psi_i,P\psi_j)\tilde{u}_1^i=\langle f,P\psi_j\rangle+\langle g_N,\psi_j\rangle_{\Gamma_N}-\langle P\tilde{g}_D,P\psi_j\rangle, \forall \psi_j\in \tilde{X}

Remembering that this formulation is for one subdomain, and in each subdomain same equations can be obtained as long as FEM is used to approximate interior solution. Take the summaion of those formulations across all the subdomains, and put them in matrix form:

\left[\begin{array}{cc}\mathbf{K_{II}} & \mathbf{K_{IC}}\\ \mathbf{K_{CI}} & \mathbf{K_{CC}} \end{array}\right]\left(\begin{array}{c}\mathbf{u_I}\\\mathbf{u_C}\end{array}\right)=\left(\begin{array}{c}\mathbf{f_I}\\\mathbf{f_C}\end{array}\right)

Which is the well-known substructure method formulation in structure engineering.  The subscripts imply “connecting” and “interior”. The unknown vectors are nodal solutions in the interior of subdomains and nodal solutions on the skeletons. Classical substructure algorithm goes like this:

  1. Find \mathbf{K_{II}}^{-1}, then \mathbf{K_{II}}^{-1}\mathbf{f_I};
  2. \mathbf{g_C}=\mathbf{f_C}-\mathbf{K_{CI}}\mathbf{K_{II}}^{-1}\mathbf{f_I};
  3. \mathbf{S_C}=\mathbf{K_C}-\mathbf{K_{CI}}\mathbf{K_{II}}^{-1}\mathbf{K_{IC}};
  4. \mathbf{u_C}=\mathbf{S_C}^{-1}\mathbf{g_C};
  5. \mathbf{u_I}=\mathbf{K_{II}}^{-1}\mathbf{f_I}-\mathbf{K_{II}}^{-1}\mathbf{K_{IC}}\mathbf{u_C}

Even though the inverse is easy to be parallelized due to the matrix being block diagonal, this algorithm is impractical for large problem since in step 4, finding schur complement is costly. Preconditioner are always used instead of solving the matrix equation in step 4.