100. FETI-DP

There is a problem with the FETI method if the bilinear-form is just \(A(u,v) = \int \nabla u \nabla v\), since then sub-domains without Dirichlet boundary conditions lead to singular \(A\)-blocks (called floating sub-domains)

The FETI - DP (where DP stands for dual/primal) is not to break all continuity, but keep some degrees of freedom continuous. E.g. keep the functions connected at the vertices. Then, the \(A\)-matrix is regular, and it is still cheap to invert.

We implement the vertex by introducing additional primal variables, and gluing the sub-domain fields to the vertex variables by a penalty method.

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())
print (mesh.GetBBoundaries())
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')
('default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default', 'default')

sub-domain spaces, plus a Vertex-space with degrees of freedom only in the sub-domain vertices. (unfortunately these freedoms are not activated per default and have to be set manually, will be fixed).

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

fesVertex = H1(mesh, definedon=mesh.BBoundaries('.*'))
fes = fes * fesVertex
print ("ndof =", fes.ndof)

freedofs = fes.FreeDofs()
for el in fes.Elements(BBND):
    freedofs.Set(el.dofs[-1])

u, v = fes.TnT()

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

ndof = 1604

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

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

uvert,vvert = u[-1], v[-1]
import ngsolve.comp
dVert = ngsolve.comp.DifferentialSymbol(BBND)

for dom in mesh.Materials('.*').Split():
    print ("vertices of dom: ", dom.Boundaries().Boundaries().Mask())
    ui, vi = domtrial[dom], domtest[dom]
    a += grad(ui)*grad(vi)*dx 
    a += 1e6*(ui-uvert)*(vi-vvert)*dVert(dom.Boundaries().Boundaries())
    f += y*x*vi*dx
    

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

vertices of dom:  0: 1111000000000000
vertices of dom:  0: 0011110000000000
vertices of dom:  0: 0000111100000000
vertices of dom:  0: 0110000011000000
vertices of dom:  0: 0010100001100000
vertices of dom:  0: 0000101000110000
vertices of dom:  0: 0000000011001100
vertices of dom:  0: 0000000001100110
vertices of dom:  0: 0000000000110011

<ngsolve.comp.LinearForm at 0x110e16670>

gfu = GridFunction(fes)

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.02465296027570176     
CG iteration 2, residual = 0.0010045080936531537     
CG iteration 3, residual = 0.0007480745192100791     
CG iteration 4, residual = 0.0006647371091368124     
CG iteration 5, residual = 0.0003811858600465746     
CG iteration 6, residual = 0.000169241279661563     
CG iteration 7, residual = 0.00017862216656462402     
CG iteration 8, residual = 0.0001139934738266756     
CG iteration 9, residual = 0.0001554495625644795     
CG iteration 10, residual = 0.000113768687735942     
CG iteration 11, residual = 2.9074188055829104e-05     
CG iteration 12, residual = 1.9459877943866364e-05     
CG iteration 13, residual = 2.6840486465394826e-05     
CG iteration 14, residual = 1.9755711867121734e-05     
CG iteration 15, residual = 1.4374652208058638e-05     
CG iteration 16, residual = 1.64519534848992e-05     
CG iteration 17, residual = 8.636931867963392e-06     
CG iteration 18, residual = 5.905016378463157e-06     
CG iteration 19, residual = 2.905220752696197e-06     
CG iteration 20, residual = 9.12863329296047e-07     
CG iteration 21, residual = 7.56247992551079e-07     
CG iteration 22, residual = 1.0041701037189113e-06     
CG iteration 23, residual = 3.1737634627058275e-07     
CG iteration 24, residual = 1.695553777750457e-07     
CG iteration 25, residual = 2.508863997641801e-07     
CG iteration 26, residual = 1.5045220307617923e-07     
CG iteration 27, residual = 7.020921520573842e-08     
CG iteration 28, residual = 4.785356420685145e-08     
CG iteration 29, residual = 3.866605392033105e-08     
CG iteration 30, residual = 2.1463017279576963e-08     
CG iteration 31, residual = 3.758459253704872e-08     
CG iteration 32, residual = 2.3586833155231812e-08     
CG iteration 33, residual = 2.5490104909186764e-08     
CG iteration 34, residual = 1.5314789582441124e-08     
CG iteration 35, residual = 5.548503269620179e-09     
CG iteration 36, residual = 3.3304662762145377e-09     
CG iteration 37, residual = 1.5999435118936547e-09     
CG iteration 38, residual = 1.9098074578346433e-09     
CG iteration 39, residual = 1.3599443467801386e-09     
CG iteration 40, residual = 8.600146982364828e-10     
CG iteration 41, residual = 9.38735155262725e-10     
CG iteration 42, residual = 5.28104194714278e-10     
CG iteration 43, residual = 3.2897454008927705e-10     
CG iteration 44, residual = 2.404058525544706e-10     
CG iteration 45, residual = 3.9906250186677225e-10     
CG iteration 46, residual = 5.337893713070697e-10     
CG iteration 47, residual = 4.0301728719558134e-10     
CG iteration 48, residual = 5.286622867929107e-10     
CG iteration 49, residual = 8.690222252824563e-10     
CG iteration 50, residual = 3.384417646519425e-09     
CG iteration 51, residual = 6.272912216339554e-09     
CG iteration 52, residual = 5.108420054998112e-09     
CG iteration 53, residual = 3.997745644108606e-09     
CG iteration 54, residual = 1.0801027340152564e-08     
CG iteration 55, residual = 8.430171598176871e-09     
CG iteration 56, residual = 8.122107471371668e-08     
CG iteration 57, residual = 5.586682956909551e-08     
CG iteration 58, residual = 6.565250197410478e-08     
CG iteration 59, residual = 9.136657458440983e-08     
CG iteration 60, residual = 8.249756181280925e-08     
CG iteration 61, residual = 1.866342862798106e-07     
CG iteration 62, residual = 1.1110095651762984e-07     
CG iteration 63, residual = 1.9794050710210138e-07     
CG iteration 64, residual = 1.54363894703238e-07     
CG iteration 65, residual = 1.2704836995365155e-07     
CG iteration 66, residual = 7.223397101208971e-08     
CG iteration 67, residual = 8.228098001464905e-08     
CG iteration 68, residual = 5.1548725281844306e-08     
CG iteration 69, residual = 2.686092528287761e-08     
CG iteration 70, residual = 8.367966401769728e-09     
CG iteration 71, residual = 7.797141147241187e-09     
CG iteration 72, residual = 5.7402344909959206e-09     
CG iteration 73, residual = 2.442091465549524e-09     
CG iteration 74, residual = 3.3710588214359143e-09     
CG iteration 75, residual = 1.4901864193034273e-09     
CG iteration 76, residual = 8.659191086951956e-10     
CG iteration 77, residual = 2.636229849355725e-10     
CG iteration 78, residual = 3.156099034312336e-10     
CG iteration 79, residual = 4.5140382881547446e-11     
CG iteration 80, residual = 1.5419782071106046e-10     
CG iteration 81, residual = 2.478080957647242e-11     
CG iteration 82, residual = 4.5891398649927797e-11     
CG iteration 83, residual = 1.59405336959108e-11     
CG iteration 84, residual = 2.945019688894114e-11     
CG iteration 85, residual = 9.804111119202031e-12     
CG iteration 86, residual = 1.2230265570371406e-11     
CG iteration 87, residual = 1.3985505538779264e-11     
CG iteration 88, residual = 4.231800553405579e-11     
CG iteration 89, residual = 2.7978007889251674e-11     
CG iteration 90, residual = 2.3626794717372264e-11     
CG iteration 91, residual = 1.1626125219562889e-10     
CG iteration 92, residual = 1.1189162193674461e-10     
CG iteration 93, residual = 5.946038351683978e-10     
CG iteration 94, residual = 9.33747991988654e-10     
CG iteration 95, residual = 1.8042842653812196e-09     
CG iteration 96, residual = 7.224245137858673e-10     
CG iteration 97, residual = 1.531974595959108e-09     
CG iteration 98, residual = 1.940887650616227e-09     
CG iteration 99, residual = 3.0438549992059032e-09     
CG iteration 100, residual = 2.231628307327755e-09     
CG iteration 101, residual = 3.313430476860122e-09     
CG iteration 102, residual = 3.3449818788895694e-09     
CG iteration 103, residual = 5.2602847475597e-09     
CG iteration 104, residual = 4.225500217695409e-09     
CG iteration 105, residual = 1.4946201253765835e-09     
CG iteration 106, residual = 1.1446064790074082e-09     
CG iteration 107, residual = 2.0174751214445433e-09     
CG iteration 108, residual = 2.7100477039474815e-09     
CG iteration 109, residual = 1.3851314670110915e-09     
CG iteration 110, residual = 3.238166250910595e-10     
CG iteration 111, residual = 2.350994115985948e-10     
CG iteration 112, residual = 1.8510294360799867e-10     
CG iteration 113, residual = 7.15901102690421e-11     
CG iteration 114, residual = 3.011034189724812e-11     
CG iteration 115, residual = 4.211033417892459e-11     
CG iteration 116, residual = 2.559910850119869e-11     
CG iteration 117, residual = 1.0099475941491787e-11     
CG iteration 118, residual = 6.357213654094172e-12     
CG iteration 119, residual = 2.003411552855566e-11     
CG iteration 120, residual = 2.0051771795038548e-12     
CG iteration 121, residual = 1.6961471486119564e-12     
CG iteration 122, residual = 1.0796321768019933e-12     
CG iteration 123, residual = 1.8268017189966915e-12     
CG iteration 124, residual = 5.010904037099035e-13     
CG iteration 125, residual = 6.238452384758804e-14     
CG iteration 126, residual = 6.667823201958212e-14     
CG iteration 127, residual = 3.78593217572049e-14     
CG iteration 128, residual = 8.617426748434672e-14     
CG iteration 129, residual = 2.0971158281186845e-14     
CG converged in 129 iterations to residual 2.0971158281186845e-14

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

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) )

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

invS = CGSolver(S, pre=pre, 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.5621483955975527     
CG iteration 2, residual = 0.025107066428877348     
CG iteration 3, residual = 0.010098600736521033     
CG iteration 4, residual = 0.0011940088258720724     
CG iteration 5, residual = 0.0003140988977799483     
CG iteration 6, residual = 4.7432636997533445e-05     
CG iteration 7, residual = 1.2451276566230542e-05     
CG iteration 8, residual = 4.135008596422449e-06     
CG iteration 9, residual = 3.1323396903690127e-06     
CG iteration 10, residual = 1.0755162319435827e-06     
CG iteration 11, residual = 2.0443270809198917e-07     
CG iteration 12, residual = 3.572652921442311e-08     
CG iteration 13, residual = 2.659552219378877e-09     
CG iteration 14, residual = 2.2763945191068638e-10     
CG iteration 15, residual = 1.4916852502521548e-11     
CG iteration 16, residual = 1.2237338394942944e-11     
CG iteration 17, residual = 6.5108492359015e-13     
CG iteration 18, residual = 4.462559659526368e-14     
CG converged in 18 iterations to residual 4.462559659526368e-14

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