98. 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 = 1612

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

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.02457114192535907     
CG iteration 2, residual = 0.0009655865792460007     
CG iteration 3, residual = 0.0007295468847148508     
CG iteration 4, residual = 0.000670213728160485     
CG iteration 5, residual = 0.0002442075276496147     
CG iteration 6, residual = 0.00018043325971218288     
CG iteration 7, residual = 0.00012931897956004144     
CG iteration 8, residual = 0.00014994098580826696     
CG iteration 9, residual = 0.0001443055151564727     
CG iteration 10, residual = 0.00011366257012602817     
CG iteration 11, residual = 3.290864111138951e-05     
CG iteration 12, residual = 1.823895122568606e-05     
CG iteration 13, residual = 2.435009942420547e-05     
CG iteration 14, residual = 2.1339462929478585e-05     
CG iteration 15, residual = 1.461279114595931e-05     
CG iteration 16, residual = 1.6048774794592202e-05     
CG iteration 17, residual = 8.915266809706737e-06     
CG iteration 18, residual = 6.142185476489926e-06     
CG iteration 19, residual = 2.896091785719707e-06     
CG iteration 20, residual = 9.464936449623437e-07     
CG iteration 21, residual = 8.284837908237672e-07     
CG iteration 22, residual = 1.0159030447596995e-06     
CG iteration 23, residual = 3.0868636483734245e-07     
CG iteration 24, residual = 1.8960711609629712e-07     
CG iteration 25, residual = 3.4765833340475694e-07     
CG iteration 26, residual = 1.6118673681228545e-07     
CG iteration 27, residual = 5.9447778758558304e-08     
CG iteration 28, residual = 6.605006365731497e-08     
CG iteration 29, residual = 4.413315859209454e-08     
CG iteration 30, residual = 2.8141302469862663e-08     
CG iteration 31, residual = 5.248474923990464e-08     
CG iteration 32, residual = 2.9043539029426944e-08     
CG iteration 33, residual = 2.2743487467415637e-08     
CG iteration 34, residual = 1.5963725052984503e-08     
CG iteration 35, residual = 4.874742605266944e-09     
CG iteration 36, residual = 2.649605943716921e-09     
CG iteration 37, residual = 2.1842565066702134e-09     
CG iteration 38, residual = 1.8185459877828952e-09     
CG iteration 39, residual = 1.6295000235476487e-09     
CG iteration 40, residual = 9.04451113893643e-10     
CG iteration 41, residual = 5.14252869167667e-10     
CG iteration 42, residual = 5.09955896632741e-10     
CG iteration 43, residual = 4.2918832827752655e-10     
CG iteration 44, residual = 2.477222913736537e-10     
CG iteration 45, residual = 4.879834751152567e-10     
CG iteration 46, residual = 5.082890761846477e-10     
CG iteration 47, residual = 3.3133477798732944e-10     
CG iteration 48, residual = 4.5113767150554525e-10     
CG iteration 49, residual = 1.095287490374656e-09     
CG iteration 50, residual = 2.145506549376324e-09     
CG iteration 51, residual = 4.262867531570349e-09     
CG iteration 52, residual = 2.3517961518027695e-09     
CG iteration 53, residual = 2.9483648024138353e-09     
CG iteration 54, residual = 8.734759226364115e-09     
CG iteration 55, residual = 7.834245127376159e-09     
CG iteration 56, residual = 6.146447866647802e-08     
CG iteration 57, residual = 4.6076691266964854e-08     
CG iteration 58, residual = 4.915824420787116e-08     
CG iteration 59, residual = 6.418534767967267e-08     
CG iteration 60, residual = 7.47962586837413e-08     
CG iteration 61, residual = 2.930553907448126e-07     
CG iteration 62, residual = 9.565798245941296e-08     
CG iteration 63, residual = 2.03636290555436e-07     
CG iteration 64, residual = 1.9712412715944527e-07     
CG iteration 65, residual = 2.530085123588813e-07     
CG iteration 66, residual = 9.272654573653059e-08     
CG iteration 67, residual = 1.499829937707831e-07     
CG iteration 68, residual = 8.120077957028037e-08     
CG iteration 69, residual = 4.3479807671432295e-08     
CG iteration 70, residual = 1.3136993255589722e-08     
CG iteration 71, residual = 1.4543206316897349e-08     
CG iteration 72, residual = 7.255924539823601e-09     
CG iteration 73, residual = 4.319768976356562e-09     
CG iteration 74, residual = 6.7291369795050536e-09     
CG iteration 75, residual = 1.5112794401643663e-09     
CG iteration 76, residual = 1.2275940107555054e-09     
CG iteration 77, residual = 1.8323475450027743e-10     
CG iteration 78, residual = 1.4065796031887282e-10     
CG iteration 79, residual = 1.0477693434382163e-10     
CG iteration 80, residual = 2.2504691139418477e-10     
CG iteration 81, residual = 4.2417977126763124e-11     
CG iteration 82, residual = 9.919608535097735e-11     
CG iteration 83, residual = 2.1062277801064125e-11     
CG iteration 84, residual = 5.344277965818667e-11     
CG iteration 85, residual = 1.7661749473608498e-11     
CG iteration 86, residual = 1.0735325656214124e-11     
CG iteration 87, residual = 1.6026463068043643e-11     
CG iteration 88, residual = 1.8214007813294235e-11     
CG iteration 89, residual = 1.1571387048872892e-11     
CG iteration 90, residual = 1.5394277134353158e-11     
CG iteration 91, residual = 4.2928674505368727e-11     
CG iteration 92, residual = 8.2316795391568e-11     
CG iteration 93, residual = 2.8959937434944356e-10     
CG iteration 94, residual = 1.1603559656307161e-09     
CG iteration 95, residual = 3.831617906213142e-10     
CG iteration 96, residual = 7.182588254146196e-10     
CG iteration 97, residual = 2.415989328156747e-09     
CG iteration 98, residual = 1.7080094873364358e-09     
CG iteration 99, residual = 1.0752970287242183e-09     
CG iteration 100, residual = 1.3192589066280015e-09     
CG iteration 101, residual = 3.787009474446208e-09     
CG iteration 102, residual = 2.94078247129712e-09     
CG iteration 103, residual = 4.9613593218524244e-09     
CG iteration 104, residual = 7.439488451851567e-09     
CG iteration 105, residual = 2.5442587575472473e-09     
CG iteration 106, residual = 3.165748555301579e-09     
CG iteration 107, residual = 4.55477547672745e-09     
CG iteration 108, residual = 3.1183020105751816e-09     
CG iteration 109, residual = 1.7708086666397952e-09     
CG iteration 110, residual = 2.4893282263672034e-10     
CG iteration 111, residual = 2.1369955072120726e-10     
CG iteration 112, residual = 3.6545462479155926e-10     
CG iteration 113, residual = 7.489646076059478e-11     
CG iteration 114, residual = 9.381250083028297e-11     
CG iteration 115, residual = 9.858021859593461e-11     
CG iteration 116, residual = 4.083426246211988e-11     
CG iteration 117, residual = 5.996530966026699e-11     
CG iteration 118, residual = 4.8010919641522714e-11     
CG iteration 119, residual = 3.2648958835637757e-11     
CG iteration 120, residual = 4.136416982885233e-12     
CG iteration 121, residual = 5.888432954887796e-12     
CG iteration 122, residual = 2.179532885997113e-12     
CG iteration 123, residual = 1.6076273998630445e-12     
CG iteration 124, residual = 2.3512726396424006e-13     
CG iteration 125, residual = 3.7726676853174435e-13     
CG iteration 126, residual = 2.330894837521124e-13     
CG iteration 127, residual = 2.0334815045743921e-13     
CG iteration 128, residual = 6.300204557356935e-14     
CG iteration 129, residual = 4.6199369110942685e-14     
CG iteration 130, residual = 1.3331751027896488e-13     
CG iteration 131, residual = 6.771162300960084e-14     
CG iteration 132, residual = 2.2207135949788015e-14     
CG converged in 132 iterations to residual 2.2207135949788015e-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.5601652811873717     
CG iteration 2, residual = 0.023767779802704447     
CG iteration 3, residual = 0.009930840510974703     
CG iteration 4, residual = 0.0011844045752136666     
CG iteration 5, residual = 0.0003179380491653574     
CG iteration 6, residual = 4.1884130318011144e-05     
CG iteration 7, residual = 1.2087414891413136e-05     
CG iteration 8, residual = 4.1781110285432425e-06     
CG iteration 9, residual = 3.6396570354424046e-06     
CG iteration 10, residual = 9.95947251901502e-07     
CG iteration 11, residual = 1.8592081188909723e-07     
CG iteration 12, residual = 2.979792841603656e-08     
CG iteration 13, residual = 2.3600669286076522e-09     
CG iteration 14, residual = 2.069381100367259e-10     
CG iteration 15, residual = 1.468752070670126e-11     
CG iteration 16, residual = 1.215306427306312e-11     
CG iteration 17, residual = 9.396284506449212e-13     
CG iteration 18, residual = 6.721280491055069e-14     
CG converged in 18 iterations to residual 6.721280491055069e-14

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