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

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.02465296027868554     
CG iteration 2, residual = 0.0010045080948046178     
CG iteration 3, residual = 0.0007480745191825382     
CG iteration 4, residual = 0.0006647371089273674     
CG iteration 5, residual = 0.00038118586162730415     
CG iteration 6, residual = 0.00016924127957501458     
CG iteration 7, residual = 0.00017862216651747526     
CG iteration 8, residual = 0.00011399347384709984     
CG iteration 9, residual = 0.00015544956276753865     
CG iteration 10, residual = 0.00011376868770978267     
CG iteration 11, residual = 2.907418805114534e-05     
CG iteration 12, residual = 1.9459877945693153e-05     
CG iteration 13, residual = 2.684048642639457e-05     
CG iteration 14, residual = 1.9755711890559478e-05     
CG iteration 15, residual = 1.4374652207083435e-05     
CG iteration 16, residual = 1.6451953487467287e-05     
CG iteration 17, residual = 8.63693187060488e-06     
CG iteration 18, residual = 5.905016379107477e-06     
CG iteration 19, residual = 2.905220753950266e-06     
CG iteration 20, residual = 9.128633295688369e-07     
CG iteration 21, residual = 7.562479949688583e-07     
CG iteration 22, residual = 1.0041704013941648e-06     
CG iteration 23, residual = 3.174160533380204e-07     
CG iteration 24, residual = 1.8539919851525282e-07     
CG iteration 25, residual = 3.1414606508220726e-07     
CG iteration 26, residual = 1.3413575487300865e-07     
CG iteration 27, residual = 7.008021756016386e-08     
CG iteration 28, residual = 4.020469339646554e-08     
CG iteration 29, residual = 4.5273708521610325e-08     
CG iteration 30, residual = 2.2452508410114673e-08     
CG iteration 31, residual = 3.703056584915259e-08     
CG iteration 32, residual = 2.2531655940300604e-08     
CG iteration 33, residual = 2.0440532749276015e-08     
CG iteration 34, residual = 1.708677864919685e-08     
CG iteration 35, residual = 3.798414063597613e-09     
CG iteration 36, residual = 4.0732143802985025e-09     
CG iteration 37, residual = 1.5162619718672763e-09     
CG iteration 38, residual = 2.2619677341594903e-09     
CG iteration 39, residual = 1.0301703999016661e-09     
CG iteration 40, residual = 9.670849241919337e-10     
CG iteration 41, residual = 6.8277295632502e-10     
CG iteration 42, residual = 6.124274268570732e-10     
CG iteration 43, residual = 3.32250192600916e-10     
CG iteration 44, residual = 2.460968479169061e-10     
CG iteration 45, residual = 6.412049294537353e-10     
CG iteration 46, residual = 3.506061968089844e-10     
CG iteration 47, residual = 4.0263310861495424e-10     
CG iteration 48, residual = 5.286614361412964e-10     
CG iteration 49, residual = 8.690230524751264e-10     
CG iteration 50, residual = 3.3807668508855315e-09     
CG iteration 51, residual = 4.478246545185634e-09     
CG iteration 52, residual = 7.786901032560145e-09     
CG iteration 53, residual = 4.467234560933242e-09     
CG iteration 54, residual = 7.62996751981704e-09     
CG iteration 55, residual = 8.159892097476285e-09     
CG iteration 56, residual = 6.310947780897848e-08     
CG iteration 57, residual = 4.205858702975159e-08     
CG iteration 58, residual = 1.2677782516426693e-07     
CG iteration 59, residual = 7.966935526980784e-08     
CG iteration 60, residual = 1.287320019479384e-07     
CG iteration 61, residual = 2.632798168202479e-07     
CG iteration 62, residual = 1.0099186473552836e-07     
CG iteration 63, residual = 2.644357965270683e-07     
CG iteration 64, residual = 1.6237942873890143e-07     
CG iteration 65, residual = 1.450884347725628e-07     
CG iteration 66, residual = 7.529716428620827e-08     
CG iteration 67, residual = 1.095681529469647e-07     
CG iteration 68, residual = 6.390790635798595e-08     
CG iteration 69, residual = 4.466917113625665e-08     
CG iteration 70, residual = 5.8496581769267444e-09     
CG iteration 71, residual = 1.0235762136016255e-08     
CG iteration 72, residual = 6.5217612495516354e-09     
CG iteration 73, residual = 2.4416372870741625e-09     
CG iteration 74, residual = 1.9983304743838603e-09     
CG iteration 75, residual = 1.8337954967115787e-09     
CG iteration 76, residual = 8.693103996879348e-10     
CG iteration 77, residual = 2.775069762006435e-10     
CG iteration 78, residual = 5.056937904632209e-10     
CG iteration 79, residual = 4.56294229882382e-11     
CG iteration 80, residual = 1.36618184731459e-10     
CG iteration 81, residual = 2.2961692695623373e-11     
CG iteration 82, residual = 4.5560709119836174e-11     
CG iteration 83, residual = 1.5256300034628712e-11     
CG iteration 84, residual = 2.084301453161241e-11     
CG iteration 85, residual = 1.0782770910841718e-11     
CG iteration 86, residual = 2.7139437932631243e-11     
CG iteration 87, residual = 9.801914203620341e-12     
CG iteration 88, residual = 5.1797649860127124e-11     
CG iteration 89, residual = 2.719447263770147e-11     
CG iteration 90, residual = 2.493416841241796e-11     
CG iteration 91, residual = 4.8707579887105684e-11     
CG iteration 92, residual = 1.084786617623542e-10     
CG iteration 93, residual = 1.3820286755443486e-09     
CG iteration 94, residual = 7.58708980075401e-10     
CG iteration 95, residual = 8.065809012485016e-10     
CG iteration 96, residual = 1.2615869374496392e-09     
CG iteration 97, residual = 2.8205626682041215e-09     
CG iteration 98, residual = 3.903338096581981e-09     
CG iteration 99, residual = 2.4463298559921924e-09     
CG iteration 100, residual = 2.3669477627155516e-09     
CG iteration 101, residual = 4.794772159980843e-09     
CG iteration 102, residual = 3.768642643037535e-09     
CG iteration 103, residual = 5.1627633120163025e-09     
CG iteration 104, residual = 4.260764059806819e-09     
CG iteration 105, residual = 2.1900967682638012e-09     
CG iteration 106, residual = 1.6889192739365799e-09     
CG iteration 107, residual = 1.1830587682732085e-09     
CG iteration 108, residual = 1.5313028199518356e-09     
CG iteration 109, residual = 3.050148731681148e-09     
CG iteration 110, residual = 4.0089900607913403e-10     
CG iteration 111, residual = 4.3052370493687505e-10     
CG iteration 112, residual = 4.1816856732960673e-10     
CG iteration 113, residual = 2.0846669478027652e-10     
CG iteration 114, residual = 3.481663622796438e-11     
CG iteration 115, residual = 2.234461405775931e-11     
CG iteration 116, residual = 4.10363491888403e-11     
CG iteration 117, residual = 8.063282306708332e-12     
CG iteration 118, residual = 7.1345178487644855e-12     
CG iteration 119, residual = 2.0181455321711325e-11     
CG iteration 120, residual = 1.993088456206386e-12     
CG iteration 121, residual = 1.7080134646151793e-12     
CG iteration 122, residual = 1.0264512411797835e-12     
CG iteration 123, residual = 1.074734390926469e-12     
CG iteration 124, residual = 7.201669157380384e-13     
CG iteration 125, residual = 2.7735933116880444e-13     
CG iteration 126, residual = 4.669430220743792e-14     
CG iteration 127, residual = 6.670758028988926e-14     
CG iteration 128, residual = 4.177453472925897e-14     
CG iteration 129, residual = 3.951870009957473e-14     
CG iteration 130, residual = 7.642934745370262e-15     
CG converged in 130 iterations to residual 7.642934745370262e-15

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.5621483956637152     
CG iteration 2, residual = 0.02510706641904706     
CG iteration 3, residual = 0.010098600733808359     
CG iteration 4, residual = 0.0011940088228786937     
CG iteration 5, residual = 0.000314098897575211     
CG iteration 6, residual = 4.7432636916910925e-05     
CG iteration 7, residual = 1.2451276569145491e-05     
CG iteration 8, residual = 4.135008604323664e-06     
CG iteration 9, residual = 3.132339685097992e-06     
CG iteration 10, residual = 1.0755162307218035e-06     
CG iteration 11, residual = 2.0443270756579144e-07     
CG iteration 12, residual = 3.5726529203105224e-08     
CG iteration 13, residual = 2.6595522191210846e-09     
CG iteration 14, residual = 2.2763945188996184e-10     
CG iteration 15, residual = 1.4916852452760398e-11     
CG iteration 16, residual = 1.2237338360609315e-11     
CG iteration 17, residual = 6.510851383162892e-13     
CG iteration 18, residual = 4.462873067049273e-14     
CG converged in 18 iterations to residual 4.462873067049273e-14

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