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 0x10b033030>

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.0840360600248355     
CG iteration 2, residual = 0.039719151026074255     
CG iteration 3, residual = 0.012296701689773551     
CG iteration 4, residual = 0.0020900138250240477     
CG iteration 5, residual = 0.0015912140518391073     
CG iteration 6, residual = 0.00031734055142703077     
CG iteration 7, residual = 0.00021581363422482137     
CG iteration 8, residual = 0.0007799867836301256     
CG iteration 9, residual = 0.0006376282665116182     
CG iteration 10, residual = 7.861701058806006e-05     
CG iteration 11, residual = 0.00010717377057099597     
CG iteration 12, residual = 4.7096192749399806e-05     
CG iteration 13, residual = 4.670865396584739e-05     
CG iteration 14, residual = 4.826747398304832e-05     
CG iteration 15, residual = 4.793177280062789e-05     
CG iteration 16, residual = 0.0003362507702467543     
CG iteration 17, residual = 3.789042259208961e-05     
CG iteration 18, residual = 0.00010052940815836863     
CG iteration 19, residual = 0.00021665830631643643     
CG iteration 20, residual = 9.238403107088956e-05     
CG iteration 21, residual = 2.7478099198131128e-05     
CG iteration 22, residual = 1.0828968506400264e-05     
CG iteration 23, residual = 3.626977130841099e-05     
CG iteration 24, residual = 1.3817656862931228e-05     
CG iteration 25, residual = 1.2579066107605714e-05     
CG iteration 26, residual = 1.2339712137642314e-05     
CG iteration 27, residual = 1.4143727527987861e-05     
CG iteration 28, residual = 8.786835855977455e-06     
CG iteration 29, residual = 1.0110142065261327e-05     
CG iteration 30, residual = 6.7128622952170054e-06     
CG iteration 31, residual = 1.9050200517461836e-05     
CG iteration 32, residual = 3.3708017175507935e-06     
CG iteration 33, residual = 6.050478593442649e-06     
CG iteration 34, residual = 2.697569661073651e-06     
CG iteration 35, residual = 8.730979542008852e-06     
CG iteration 36, residual = 2.2746981095664904e-06     
CG iteration 37, residual = 6.613459255455237e-06     
CG iteration 38, residual = 1.0748016619442559e-06     
CG iteration 39, residual = 5.885165437134679e-07     
CG iteration 40, residual = 6.747146301790482e-07     
CG iteration 41, residual = 2.088423846772488e-06     
CG iteration 42, residual = 4.966642520897391e-07     
CG iteration 43, residual = 4.934588420886047e-07     
CG iteration 44, residual = 8.550249913595322e-07     
CG iteration 45, residual = 4.6680467729668795e-07     
CG iteration 46, residual = 3.634358170841065e-07     
CG iteration 47, residual = 5.032769912102492e-07     
CG iteration 48, residual = 2.830769329908982e-07     
CG iteration 49, residual = 7.042066167756622e-07     
CG iteration 50, residual = 1.501364906305562e-06     
CG iteration 51, residual = 3.3846079792464547e-07     
CG iteration 52, residual = 2.589184354351973e-07     
CG iteration 53, residual = 2.2425254228350306e-07     
CG iteration 54, residual = 7.42925292753348e-07     
CG iteration 55, residual = 9.083985320499209e-08     
CG iteration 56, residual = 6.780957388869106e-08     
CG iteration 57, residual = 5.770415667746401e-07     
CG iteration 58, residual = 3.242815265591493e-08     
CG iteration 59, residual = 2.699649464325231e-08     
CG iteration 60, residual = 2.4993800362231312e-08     
CG iteration 61, residual = 1.7370248588086147e-07     
CG iteration 62, residual = 3.280837377997276e-08     
CG iteration 63, residual = 3.190275311155329e-08     
CG iteration 64, residual = 3.409472839829506e-08     
CG iteration 65, residual = 1.807372753490529e-08     
CG iteration 66, residual = 1.6843054543068204e-08     
CG iteration 67, residual = 4.454894178581335e-08     
CG iteration 68, residual = 6.477665274484601e-09     
CG iteration 69, residual = 9.980680775081782e-09     
CG iteration 70, residual = 6.57716816137978e-09     
CG iteration 71, residual = 6.3818103489254606e-09     
CG iteration 72, residual = 5.8810371594804925e-09     
CG iteration 73, residual = 4.943124857589382e-09     
CG iteration 74, residual = 1.2756432221338466e-08     
CG iteration 75, residual = 4.149090991030669e-09     
CG iteration 76, residual = 7.820778139584721e-09     
CG iteration 77, residual = 1.8576707725878566e-09     
CG iteration 78, residual = 3.2007754421720098e-09     
CG iteration 79, residual = 1.51062684534242e-09     
CG iteration 80, residual = 1.2294631506245635e-09     
CG iteration 81, residual = 2.3736741231830405e-09     
CG iteration 82, residual = 9.905233480851317e-10     
CG iteration 83, residual = 1.5246855534366472e-09     
CG iteration 84, residual = 9.805764468819512e-10     
CG iteration 85, residual = 4.909571309816757e-09     
CG iteration 86, residual = 1.4594671894532613e-09     
CG iteration 87, residual = 8.779053607170801e-10     
CG iteration 88, residual = 2.566824448448467e-10     
CG iteration 89, residual = 1.8722139367717476e-09     
CG iteration 90, residual = 2.0767828255207989e-10     
CG iteration 91, residual = 1.3199099234624887e-10     
CG iteration 92, residual = 3.737762641254479e-10     
CG iteration 93, residual = 7.879596270824799e-11     
CG iteration 94, residual = 7.088299692363708e-10     
CG iteration 95, residual = 2.327669068256477e-10     
CG iteration 96, residual = 1.5105967773065588e-10     
CG iteration 97, residual = 1.2225551315020386e-10     
CG iteration 98, residual = 5.1102696147400337e-11     
CG iteration 99, residual = 3.223218202034396e-11     
CG iteration 100, residual = 4.589681348912029e-11     
CG iteration 101, residual = 2.9176979372452557e-11     
CG iteration 102, residual = 2.7087303748416527e-11     
CG iteration 103, residual = 3.3014230754789665e-11     
CG iteration 104, residual = 8.296943858583755e-12     
CG iteration 105, residual = 4.17629892856994e-11     
CG iteration 106, residual = 4.157455547812243e-11     
CG iteration 107, residual = 1.1743741659880897e-11     
CG iteration 108, residual = 2.759853237700843e-11     
CG iteration 109, residual = 1.5973677258816093e-11     
CG iteration 110, residual = 3.6419178667718255e-12     
CG iteration 111, residual = 5.271750646227386e-12     
CG iteration 112, residual = 6.6272890569287374e-12     
CG iteration 113, residual = 3.6140703231896918e-12     
CG iteration 114, residual = 2.9915748340383654e-12     
CG iteration 115, residual = 1.3043769318572222e-11     
CG iteration 116, residual = 3.409623094499772e-12     
CG iteration 117, residual = 1.6501859226815918e-12     
CG iteration 118, residual = 1.3639582257123893e-12     
CG iteration 119, residual = 2.0101001437153654e-12     
CG iteration 120, residual = 7.503944067158417e-13     
CG iteration 121, residual = 2.14590222827916e-12     
CG iteration 122, residual = 9.04926769077577e-13     
CG iteration 123, residual = 3.972178431037747e-13     
CG iteration 124, residual = 6.490950940610273e-13     
CG iteration 125, residual = 5.00914946858202e-13     
CG iteration 126, residual = 1.962022711980884e-12     
CG iteration 127, residual = 3.6292936604710116e-13     
CG iteration 128, residual = 7.171620363509996e-13     
CG iteration 129, residual = 4.101543108685856e-13     
CG iteration 130, residual = 5.817695605321162e-13     
CG iteration 131, residual = 2.7707678391773793e-13     
CG iteration 132, residual = 1.9899200352689228e-13     
CG iteration 133, residual = 9.444949943002613e-13     
CG iteration 134, residual = 1.5013435065561195e-13     
CG iteration 135, residual = 4.0642647152259627e-13     
CG iteration 136, residual = 2.400963777161998e-13     
CG iteration 137, residual = 4.489639015937824e-14     
CG converged in 137 iterations to residual 4.489639015937824e-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/2dc70abe4d033d57815f3624aaa5b8c184ed34feac2cedb580124daacad88297.png

<Figure size 640x480 with 0 Axes>

Draw(gftot, mesh);

FETI methods

Contents

99. FETI methods#

99.1. Preconditioner for \(S\)#