FETI methods

99. FETI methods#

Finite Element Tearing and Interconnection.

As the name says, we break the global (finite element) system apart, and then reenforce continuity via Lagrange parameters.

from ngsolve import *
from ngsolve.webgui import Draw

from netgen.occ import *
shape = None
mx, my = 3,3
for i in range(mx): 
    for j in range(my):
        rect = MoveTo(i/mx, j/mx).Rectangle(1/mx, 1/my).Face()
        rect.faces.name='mat'+str(i)+str(j)
        shape = Glue([shape,rect]) if shape else rect

shape.edges[Y<1e-4].name = "bot"
mesh = Mesh(OCCGeometry(shape, dim=2).GenerateMesh(maxh=0.04))
print (mesh.GetMaterials())
print (mesh.GetBoundaries())
Draw (mesh);

('mat00', 'mat01', 'mat02', 'mat10', 'mat11', 'mat12', 'mat20', 'mat21', 'mat22')
('bot', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'bot', 'default', 'default', 'default', 'default', 'default', 'default', 'bot', 'default', 'default', 'default', 'default', 'default', 'default')

We define \(H^1\)-spaces on sub-domains \(\Omega_i\). The product spaces is

\[ V = \Pi_i H^1(\Omega_i) \]

fes = None
for dom in mesh.Materials('.*').Split():
    fesi = Compress(H1(mesh, definedon=dom, dirichlet="bot"))
    fes = fes * fesi if fes else fesi

print ("ndof =", fes.ndof)

u, v = fes.TnT()

domtrial = {} 
domtest = {}
for nr,dom in enumerate (mesh.Materials('.*').Split()):
    domtrial[dom] = u[nr]
    domtest[dom] = v[nr]

ndof = 854

We identify the interfaces (i.e. internal boundaries) between two sub-domains. On these \(\gamma_{ij}\) we define spaces for the Lagrange parameter. Although the spaces are \(H^{-1/2}\), we use \(H^1(\gamma_{ij})\) to obtain the same number of constraints as we have basis functions on the interface.

feslam = None
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        feslami = Compress(H1(mesh, definedon=inter))
        feslam = feslam * feslami if feslam else feslami 
    
print ("ndof-lam:", feslam.ndof)

lam, mu = feslam.TnT()

intertrial = {} 
intertest = {}
nr = 0
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        intertrial[inter] = lam[nr]
        intertest[inter] = mu[nr]
        nr = nr+1

ndof-lam: 108

We define bilinear-forms on the sub-domains:

\[ a(u,v) = \sum_{\Omega_i} \int_{\Omega_i} \nabla u_i \nabla v_i + u_i v_i \]

and the constraint equations

\[ b(u,\mu) = \sum_{\gamma_{ij}} \int_{\gamma_{ij}} (u_i - u_j) \mu_{ij} \]

a = BilinearForm(fes)
f = LinearForm(fes)
b = BilinearForm(trialspace=fes, testspace=feslam)

for (ui,vi) in zip(u,v):
    a += grad(ui)*grad(vi)*dx + 1*ui*vi*dx
    f += y*x*vi*dx
    

for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        dom1,dom2 = doms
        # a += 1*(domtrial[dom1]-domtrial[dom2])*(domtest[dom1]-domtest[dom2])*ds(inter)
        b += (domtrial[dom1]-domtrial[dom2]) * intertest[inter] * ds(inter)
        
a.Assemble()
b.Assemble()
f.Assemble()

<ngsolve.comp.LinearForm at 0x10c72e1b0>

Obviously, if we only solve the decomposed sub-problems, we don’t get the correct solution:

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
gftot = CF ( list(gfu.components) )
Draw (gftot, mesh);

Next we solve the saddle-point problem

\[\begin{split} \left( \begin{array}{cc} A & B^T \\ B & 0 \end{array} \right) \left( \begin{array}{c} u \\ \lambda \end{array} \right) = \left( \begin{array}{c} f \\ 0 \end{array} \right) \end{split}\]

We explicitly build the Schur-complement \(S = B A^{-1} B^T\), and use conjugate gradients to solve for the Lagrange parameter \(\lambda\). Then we recover the primal variable \(u\) from the first equation:

ainv = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())
S = b.mat @ ainv @ b.mat.T
g = (b.mat @ ainv * f.vec).Evaluate()

from ngsolve.krylovspace import CGSolver
invS = CGSolver(S, pre=IdentityMatrix(feslam.ndof), printrates="\r", maxiter=500)

lam = g.CreateVector()
lam.data = invS * g
gfu.vec.data = ainv * (f.vec - b.mat.T * lam)

CG iteration 1, residual = 0.08403606002482938     
CG iteration 2, residual = 0.03971915102607125     
CG iteration 3, residual = 0.012296701689717297     
CG iteration 4, residual = 0.0020900138250192013     
CG iteration 5, residual = 0.001591214051841302     
CG iteration 6, residual = 0.0003173405514275741     
CG iteration 7, residual = 0.0002158136440848936     
CG iteration 8, residual = 0.000779986340623723     
CG iteration 9, residual = 0.0006376282700496858     
CG iteration 10, residual = 7.861701055840406e-05     
CG iteration 11, residual = 0.00010717377057098683     
CG iteration 12, residual = 4.709619274837254e-05     
CG iteration 13, residual = 4.670865396560278e-05     
CG iteration 14, residual = 4.826747522390014e-05     
CG iteration 15, residual = 4.793910703776113e-05     
CG iteration 16, residual = 0.00033432359906604653     
CG iteration 17, residual = 4.040105453373295e-05     
CG iteration 18, residual = 7.761223966339449e-05     
CG iteration 19, residual = 0.00017044087576055045     
CG iteration 20, residual = 0.0001823362812443871     
CG iteration 21, residual = 2.6581938437272876e-05     
CG iteration 22, residual = 1.0468586135634859e-05     
CG iteration 23, residual = 7.075179285880117e-05     
CG iteration 24, residual = 9.851534588915696e-06     
CG iteration 25, residual = 2.0337433831242728e-05     
CG iteration 26, residual = 1.364602798153885e-05     
CG iteration 27, residual = 1.4148871295746056e-05     
CG iteration 28, residual = 8.78788708389717e-06     
CG iteration 29, residual = 2.4681501328341185e-05     
CG iteration 30, residual = 7.149810128775211e-06     
CG iteration 31, residual = 8.617938734383813e-06     
CG iteration 32, residual = 3.5961200258358228e-06     
CG iteration 33, residual = 9.194640649505682e-06     
CG iteration 34, residual = 3.147205177517952e-06     
CG iteration 35, residual = 3.6508823019970142e-06     
CG iteration 36, residual = 2.1714062809925087e-06     
CG iteration 37, residual = 4.5749957296047e-06     
CG iteration 38, residual = 1.1024937934305609e-06     
CG iteration 39, residual = 5.878799317838693e-07     
CG iteration 40, residual = 6.782556290526538e-07     
CG iteration 41, residual = 2.2756154441538703e-06     
CG iteration 42, residual = 4.965205503262072e-07     
CG iteration 43, residual = 5.546959504584651e-07     
CG iteration 44, residual = 5.837520698945262e-07     
CG iteration 45, residual = 4.57845697728726e-07     
CG iteration 46, residual = 7.231835673812483e-07     
CG iteration 47, residual = 5.439503474479101e-07     
CG iteration 48, residual = 4.320265453489732e-07     
CG iteration 49, residual = 2.648855387948286e-07     
CG iteration 50, residual = 1.134425614201421e-06     
CG iteration 51, residual = 2.413224961647785e-07     
CG iteration 52, residual = 3.949443213808673e-07     
CG iteration 53, residual = 5.267482364906711e-07     
CG iteration 54, residual = 2.3511076761996342e-07     
CG iteration 55, residual = 9.083818889508835e-08     
CG iteration 56, residual = 6.773574769704733e-08     
CG iteration 57, residual = 3.779089351699628e-07     
CG iteration 58, residual = 3.2504798645334696e-08     
CG iteration 59, residual = 3.181058264569693e-08     
CG iteration 60, residual = 5.0183694921978756e-08     
CG iteration 61, residual = 2.4964214374600503e-08     
CG iteration 62, residual = 4.726478357917658e-08     
CG iteration 63, residual = 2.5039047312674654e-08     
CG iteration 64, residual = 3.845483865103992e-08     
CG iteration 65, residual = 1.3559047885475959e-08     
CG iteration 66, residual = 2.7214372521195454e-08     
CG iteration 67, residual = 8.584129886287469e-09     
CG iteration 68, residual = 9.681202723842704e-09     
CG iteration 69, residual = 5.058481492802685e-09     
CG iteration 70, residual = 3.173029753095916e-08     
CG iteration 71, residual = 1.022984964488607e-08     
CG iteration 72, residual = 4.275550523956888e-09     
CG iteration 73, residual = 1.2097122378858898e-08     
CG iteration 74, residual = 7.2612754947924205e-09     
CG iteration 75, residual = 1.243275716288881e-08     
CG iteration 76, residual = 1.7899947787476044e-09     
CG iteration 77, residual = 2.3282639731436342e-09     
CG iteration 78, residual = 1.743160715384927e-09     
CG iteration 79, residual = 3.710268657784957e-09     
CG iteration 80, residual = 4.377257749189571e-09     
CG iteration 81, residual = 1.2736523027834313e-09     
CG iteration 82, residual = 7.879091085908694e-10     
CG iteration 83, residual = 3.98516755093122e-09     
CG iteration 84, residual = 9.225889321774107e-10     
CG iteration 85, residual = 1.4987487635092555e-09     
CG iteration 86, residual = 9.720520337689996e-10     
CG iteration 87, residual = 2.936306052014446e-10     
CG iteration 88, residual = 5.664529166492434e-10     
CG iteration 89, residual = 5.531303118707229e-10     
CG iteration 90, residual = 2.1913571764242577e-10     
CG iteration 91, residual = 1.248919249845544e-10     
CG iteration 92, residual = 9.386237196235745e-10     
CG iteration 93, residual = 7.808585810456647e-11     
CG iteration 94, residual = 4.657764956730074e-10     
CG iteration 95, residual = 1.38423407946859e-10     
CG iteration 96, residual = 2.514193719549919e-10     
CG iteration 97, residual = 5.364089874199844e-10     
CG iteration 98, residual = 4.805630947241368e-11     
CG iteration 99, residual = 5.103354589997376e-11     
CG iteration 100, residual = 3.0815055730829476e-11     
CG iteration 101, residual = 1.2324868434090698e-10     
CG iteration 102, residual = 2.0099103452177086e-11     
CG iteration 103, residual = 9.368650549763222e-12     
CG iteration 104, residual = 8.957734864111044e-12     
CG iteration 105, residual = 3.1650803276447296e-11     
CG iteration 106, residual = 5.107488947912795e-11     
CG iteration 107, residual = 6.108746630352481e-12     
CG iteration 108, residual = 9.793646974287008e-12     
CG iteration 109, residual = 2.218510650706047e-11     
CG iteration 110, residual = 5.155332871951183e-12     
CG iteration 111, residual = 6.664391434870093e-12     
CG iteration 112, residual = 4.485631530525624e-12     
CG iteration 113, residual = 2.601635370497023e-12     
CG iteration 114, residual = 5.368450584743921e-12     
CG iteration 115, residual = 1.745779868975074e-11     
CG iteration 116, residual = 1.7548933669526746e-12     
CG iteration 117, residual = 1.5294775136429558e-12     
CG iteration 118, residual = 2.083312040727867e-12     
CG iteration 119, residual = 6.926357287435466e-13     
CG iteration 120, residual = 2.1455178979036545e-12     
CG iteration 121, residual = 9.499168518365355e-13     
CG iteration 122, residual = 4.146481710994661e-12     
CG iteration 123, residual = 3.4520879930032936e-13     
CG iteration 124, residual = 1.4004478124175422e-12     
CG iteration 125, residual = 4.710136319442546e-13     
CG iteration 126, residual = 3.3697220676569303e-13     
CG iteration 127, residual = 2.9845076689035364e-13     
CG iteration 128, residual = 7.049554175392313e-13     
CG iteration 129, residual = 3.9912445057541187e-13     
CG iteration 130, residual = 8.014346198864077e-13     
CG iteration 131, residual = 5.602405682920977e-13     
CG iteration 132, residual = 2.1806352801243222e-13     
CG iteration 133, residual = 1.449388145263939e-13     
CG iteration 134, residual = 3.993028646781513e-13     
CG iteration 135, residual = 6.000634588011272e-13     
CG iteration 136, residual = 2.1616313787824553e-13     
CG iteration 137, residual = 4.3463448709645174e-14     
CG converged in 137 iterations to residual 4.3463448709645174e-14

gftot = CF ( list(gfu.components) )
Draw(gftot, mesh);

99.1. Preconditioner for \(S\)#

The function space for the Lagrange parameter is \(\Pi_{ij} H^{-1/2} (\gamma_{ij})\). The Schur complement matrix scales like a bilinear-form in \(H^{-1/2}\). As a preconditioner, we need a matrix which scales like \(\Pi_{ij} H^{1/2}\). We cheat a bit with the non-additivity of the \(H^{1/2}\)-norm, and use \(\sum_{\Omega_i} \| \operatorname{tr} u \|_{H^{1/2}(\partial \Omega_i)}^2\)

We have to map functions from the skeleton onto the domain:

bnddofs = fes.GetDofs(mesh.Boundaries(".*"))
innerdofs = ~bnddofs

massbnd = BilinearForm(fes)
for (ui,vi) in zip(u,v):
    massbnd += ui*vi*ds
massbnd.Assemble()
invmassbnd = massbnd.mat.Inverse(inverse="sparsecholesky", freedofs=bnddofs)

massinter = BilinearForm(feslam)
for inter in mesh.Boundaries('.*').Split():
    doms = inter.Neighbours(VOL).Split()
    if len(doms) == 2:
        massinter += intertrial[inter]*intertest[inter]*ds(inter)
massinter.Assemble()

emb = invmassbnd@b.mat.T@massinter.mat.Inverse(inverse="sparsecholesky")

used dof inconsistency
(silence this warning by setting BilinearForm(...check_unused=False) )

The \(H^{1/2}(\partial \Omega_i)\)-norms are obtained by forming Schur-complements of the sub-domain matrices with respect to boundary dofs:

SchurDir = a.mat - a.mat@a.mat.Inverse(inverse="sparsecholesky", freedofs=innerdofs)@a.mat

pre = emb.T @ SchurDir @ emb

invS = solvers.CGSolver(S, pre=pre, plotrates=True, maxiter=500)

lam = g.CreateVector()
lam.data = invS * g
gfu.vec.data = ainv * (f.vec - b.mat.T * lam)

../_images/25ff51885c5c80891dbc705dfcf970e62ae963d4425f7a2ba8527fc47721b75a.png

<Figure size 640x480 with 0 Axes>

Draw(gftot, mesh);

FETI methods

Contents

99. FETI methods#

99.1. Preconditioner for \(S\)#