# Composing a solution from subregions

Hello everyone; I would like to be able to split my simulation region into parts; solve those partially, combine those to get the solution.
The reason is that the problem is composed of repeating structures and I hope to increase the calculation speed.
Can anyone point to an example; I could not find any.

Have a look into this:

https://jschoeberl.github.io/iFEM/subspacecorrection/minimaldd.html

and also section 72 ff

Hope this gives some ideas.

Joachim

Thank you very much; I was not aware of that resource. That should at least keep me busy for some time.

I had a look at the given example. Since I am using a more complex geometry I tried to extend the example. As a first step I wanted to combine sets of 3 rectangles to form 3 regions.

``````    regions = []
key = lambda x: x[-1]
for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
reg = reduce(operator.add, (mesh.Materials(m) for m in g))
regions.append(reg)

fes = FESpace([Compress(H1(mesh, definedon=dom, dirichlet="bot")) for dom in regions])

a = BilinearForm(fes)
f = LinearForm(fes)

for (ui,vi) in zip(u,v):
f += y*x*vi*dx

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
gftot = CF(list(gfu.components))
``````

This clearly does not work; (simply check by actually composing the the geometry to 3 larger rectangles).
What am I missing here.
Best to everyone

if you post a complete example, your chance for help increases
Joachim

My apologies; posting a minimal working example should have been a no-brainer.

``````
from ngsolve import *
from netgen.geom2d import CSG2d, Rectangle
from itertools import groupby
from functools import reduce
import operator

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)/my), \
mat='mat'+str(i)+str(j), \
bc = 'default', \
bottom = 'bot' if j == 0 else None)

mesh = Mesh(geo.GenerateMesh(maxh=0.04))

key = lambda x: x[-1]

mats = ['|'.join(m for m in g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]
regions = [reg:=mesh.Materials(mat) for mat in mats]

# equivalent result using this construction of the regions
# mats = [list(g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]
# regions = []
# for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
#     reg = reduce(operator.add, (mesh.Materials(m) for m in g))
#     regions.append(reg)

fes = FESpace([Compress(H1(mesh, definedon=dom, dirichlet="bot")) for dom in regions])

print ("ndof =", fes.ndof)

u, v = fes.TnT()

a = BilinearForm(fes)
f = LinearForm(fes)

for (ui,vi) in zip(u,v):
f += y*x*vi*dx

a.Assemble()
f.Assemble()

gfu = GridFunction(fes)
gfu.vec.data = a.mat.Inverse(inverse="sparsecholesky", freedofs=fes.FreeDofs())*f.vec
gftot = CF(list(gfu.components))
Draw(gftot)
``````

Using the method mentioned in your initial response, I can get it to work.
Still, I am confused what went wrong with the other approach

``````
from ngsolve import *
from netgen.geom2d import CSG2d, Rectangle
from itertools import groupby
from functools import reduce
import operator
from ngsolve.krylovspace import CGSolver

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)/my), \
mat='mat'+str(i)+str(j), \
bc = 'default', \
bottom = 'bot' if j == 0 else None)

mesh = Mesh(geo.GenerateMesh(maxh=0.04))

key = lambda x: x[-1]
mats = [list(g) for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key)]

regions = []
for k, g in groupby(sorted(mesh.GetMaterials(), key=key), key):
reg = reduce(operator.add, (mesh.Materials(m) for m in g))
regions.append(reg)

fes = H1(mesh, order=1, dirichlet='bot')

print ("ndof =", fes.ndof)

u, v = fes.TnT()

a = BilinearForm(fes)
f = LinearForm(fes)
f += x*v*dx
f.Assemble()
a.Assemble()

invs = []
for dom in regions:
domaindofs = fes.GetDofs(dom) & fes.FreeDofs()
print ("num dofs:", domaindofs.NumSet())
invs.append(a.mat.Inverse(domaindofs))

inv = CGSolver(mat=a.mat, pre=invs)
print ("Iterations: ", len(inv.errors))
gfu = GridFunction(fes)
gfu.vec.data = inv * f.vec
gftot = CF(gfu)
Draw(gftot)
``````

Some apologies at this point: I got it to work; since I use custom drawing functions, drawing the result was wrong. it is still a bit surprising but with the following replacement I got it to work:

``````gftot = reduce(operator.add, (CF(comp) for comp in gfu.components))