81. 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.geom2d import CSG2d, Rectangle
geo = CSG2d()
mx, my = 3,3
for i in range(mx):
for j in range(my):
rect = Rectangle(pmin=(i/mx,j/my), \
pmax=((i+1)/mx,(j+1)/mx), \
mat='mat'+str(i)+str(j), \
bc = 'default', \
bottom = 'bot' if j == 0 else None)
geo.Add(rect)
mesh = Mesh(geo.GenerateMesh(maxh=0.05))
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', 'default', 'default', 'bot', 'default', 'default', 'default', 'default', 'default', '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', '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 = 1110
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: 96
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: 0011010100000000
vertices of dom: 0: 0000010110100000
vertices of dom: 0: 0101101000000000
vertices of dom: 0: 0001011001000000
vertices of dom: 0: 0000010011010000
vertices of dom: 0: 0000101000001100
vertices of dom: 0: 0000001001000110
vertices of dom: 0: 0000000001010011
<ngsolve.comp.LinearForm at 0x10e19adf0>
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.02435467077440151
CG iteration 2, residual = 0.0010069624297562696
CG iteration 3, residual = 0.000725752752707929
CG iteration 4, residual = 0.0007372917166384283
CG iteration 5, residual = 0.00048575900897660407
CG iteration 6, residual = 0.00013426732432066454
CG iteration 7, residual = 0.00013119349751905874
CG iteration 8, residual = 9.75197995777629e-05
CG iteration 9, residual = 0.00015055318197175534
CG iteration 10, residual = 7.97741566142421e-05
CG iteration 11, residual = 2.5404009954080077e-05
CG iteration 12, residual = 1.132574138753976e-05
CG iteration 13, residual = 1.0260991546942733e-05
CG iteration 14, residual = 5.181919481874446e-06
CG iteration 15, residual = 2.13029395438108e-06
CG iteration 16, residual = 1.7112644172491408e-06
CG iteration 17, residual = 8.984965667907077e-07
CG iteration 18, residual = 4.938071598197655e-07
CG iteration 19, residual = 6.285198733015978e-07
CG iteration 20, residual = 4.687990675892528e-07
CG iteration 21, residual = 1.6687818393381375e-07
CG iteration 22, residual = 1.0305053745835908e-07
CG iteration 23, residual = 7.078183206613188e-08
CG iteration 24, residual = 5.404563477407528e-08
CG iteration 25, residual = 6.692753635219975e-08
CG iteration 26, residual = 2.1237307982241775e-08
CG iteration 27, residual = 1.2924811429051319e-08
CG iteration 28, residual = 6.924816628679435e-09
CG iteration 29, residual = 3.858427395202262e-09
CG iteration 30, residual = 1.738830450395595e-09
CG iteration 31, residual = 3.269058169182027e-09
CG iteration 32, residual = 1.4731343120998787e-09
CG iteration 33, residual = 3.2755319017225014e-09
CG iteration 34, residual = 1.1363485282619605e-09
CG iteration 35, residual = 1.055292332179194e-09
CG iteration 36, residual = 4.373443200688844e-10
CG iteration 37, residual = 2.8470072131627334e-10
CG iteration 38, residual = 4.746744809618211e-10
CG iteration 39, residual = 4.848906170432547e-10
CG iteration 40, residual = 1.1080931524797886e-09
CG iteration 41, residual = 9.021680172914178e-10
CG iteration 42, residual = 1.714008685346325e-09
CG iteration 43, residual = 1.843210063311627e-09
CG iteration 44, residual = 2.809635276715441e-09
CG iteration 45, residual = 9.489842622889638e-09
CG iteration 46, residual = 3.3357651872399246e-08
CG iteration 47, residual = 2.1690509397482924e-08
CG iteration 48, residual = 1.859291633883496e-08
CG iteration 49, residual = 3.343495008483151e-08
CG iteration 50, residual = 7.697808726680315e-08
CG iteration 51, residual = 2.1168648242240292e-07
CG iteration 52, residual = 1.626339847741722e-07
CG iteration 53, residual = 3.5530923154043247e-07
CG iteration 54, residual = 6.207797463801602e-08
CG iteration 55, residual = 8.506336456463105e-08
CG iteration 56, residual = 5.096870591581053e-08
CG iteration 57, residual = 3.285261803064561e-08
CG iteration 58, residual = 5.3782375941630555e-08
CG iteration 59, residual = 1.0938448611287085e-08
CG iteration 60, residual = 1.5142198151406168e-08
CG iteration 61, residual = 1.6628845003728513e-08
CG iteration 62, residual = 3.4639916534842686e-09
CG iteration 63, residual = 5.3017815971438115e-09
CG iteration 64, residual = 6.533828877556726e-09
CG iteration 65, residual = 2.8322959488698175e-09
CG iteration 66, residual = 4.49950661134314e-09
CG iteration 67, residual = 2.345286858268262e-10
CG iteration 68, residual = 3.7584634540051176e-10
CG iteration 69, residual = 7.268337563450049e-11
CG iteration 70, residual = 9.839719533898004e-11
CG iteration 71, residual = 4.157232828249118e-11
CG iteration 72, residual = 1.6852922515902963e-11
CG iteration 73, residual = 1.557134991353894e-11
CG iteration 74, residual = 2.0401019232202664e-11
CG iteration 75, residual = 2.5515751696702515e-11
CG iteration 76, residual = 2.326469974823354e-11
CG iteration 77, residual = 1.032001079597646e-10
CG iteration 78, residual = 1.0197734469007467e-10
CG iteration 79, residual = 6.627318460818314e-10
CG iteration 80, residual = 3.391198049839224e-10
CG iteration 81, residual = 2.473282009251378e-09
CG iteration 82, residual = 7.806223719927285e-10
CG iteration 83, residual = 7.0871232244493624e-09
CG iteration 84, residual = 6.1750190919030915e-09
CG iteration 85, residual = 6.967816862422729e-09
CG iteration 86, residual = 7.879557505654338e-09
CG iteration 87, residual = 2.043410523067387e-09
CG iteration 88, residual = 3.586717823012357e-09
CG iteration 89, residual = 3.389007726660379e-09
CG iteration 90, residual = 2.7766245457196026e-09
CG iteration 91, residual = 9.273342095791625e-10
CG iteration 92, residual = 3.9597946922305316e-10
CG iteration 93, residual = 2.3106696194022617e-10
CG iteration 94, residual = 4.5767269529447143e-10
CG iteration 95, residual = 1.6886868114230752e-10
CG iteration 96, residual = 7.808004503566659e-11
CG iteration 97, residual = 1.097469463440972e-10
CG iteration 98, residual = 6.984341809190383e-11
CG iteration 99, residual = 2.0861189697629896e-11
CG iteration 100, residual = 7.56906916257152e-12
CG iteration 101, residual = 4.0349543337339e-12
CG iteration 102, residual = 1.854864477206059e-12
CG iteration 103, residual = 3.803880010352867e-12
CG iteration 104, residual = 1.2322662719717318e-12
CG iteration 105, residual = 1.0665306914429312e-12
CG iteration 106, residual = 1.021302779591812e-12
CG iteration 107, residual = 4.2778255546402775e-13
CG iteration 108, residual = 1.6904904901851865e-13
CG iteration 109, residual = 8.837272655973157e-14
CG iteration 110, residual = 2.5871503487735405e-14
CG iteration 111, residual = 1.2079866858818379e-14
CG converged in 111 iterations to residual 1.2079866858818379e-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.5225189901845548
CG iteration 2, residual = 0.01919012190952205
CG iteration 3, residual = 0.007175095266953569
CG iteration 4, residual = 0.0007041293195152594
CG iteration 5, residual = 0.00014885168670214898
CG iteration 6, residual = 2.7978921512925343e-05
CG iteration 7, residual = 7.323436571077358e-06
CG iteration 8, residual = 5.782735872932991e-06
CG iteration 9, residual = 4.253381714290263e-06
CG iteration 10, residual = 6.042700555926476e-07
CG iteration 11, residual = 1.8920270278170713e-07
CG iteration 12, residual = 1.1315250560178633e-08
CG iteration 13, residual = 8.756363093945918e-10
CG iteration 14, residual = 8.434177851031543e-11
CG iteration 15, residual = 5.834859801621774e-11
CG iteration 16, residual = 5.752780375295666e-12
CG iteration 17, residual = 1.7683333378414936e-13
CG converged in 17 iterations to residual 1.7683333378414936e-13
gftot = CF ( list(gfu.components) )
Draw(gftot, mesh);