This page was generated from unit-1.5-subdomains/subdomains.ipynb.

1.5 Spaces and forms on subdomains¶

In NGSolve, finite element spaces can be defined on subdomains. This is useful for multiphysics problems like fluid-structure interaction.

In addition, bilinear or linear forms can be defined as integrals over regions. Regions are parts of the domain. They may be subdomains, or parts of the domain boundary, or parts of subdomain interfaces.

In this tutorial, you will learn about

defining finite element spaces on Regions,
Compressing such spaces,
integrating on Regionss, and
set operations on Region objects via Mask.

[1]:

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *   # Opencascade for geometry modeling

Naming subdomains and boundaries¶

We define a geometry with multiple regions and assign names to these regions below:

[2]:

outer = Rectangle(2, 2).Face()
outer.edges.name="outer"
outer.edges.Max(X).name = "r"
outer.edges.Min(X).name = "l"
outer.edges.Min(Y).name = "b"
outer.edges.Max(Y).name = "t"

inner = MoveTo(1, 0.9).Rectangle(0.3, 0.5).Face()
inner.edges.name="interface"
outer = outer - inner

inner.faces.name="inner"
inner.faces.col = (1, 0, 0)
outer.faces.name="outer"

geo = Glue([inner, outer])
Draw(geo);

[3]:

mesh = Mesh(OCCGeometry(geo, dim=2).GenerateMesh(maxh=0.2))
Draw(mesh);

This mesh has two subdomains that we named "inner" and "outer" during the geometry construction. They can be identified as Region objects.

[4]:

mesh.Materials("inner")

[4]:

<ngsolve.comp.Region at 0x7fec44ed5db0>

The bottom, right, top and left parts of the outer rectangle’s boundaries define boundary segments, which we respectively labeled "b", "r", "t", "l". They are also instances of Region objects:

[5]:

mesh.Boundaries("b")

[5]:

<ngsolve.comp.Region at 0x7fec44ed55b0>

A finite element space on a subdomain¶

[6]:

fes1 = H1(mesh, definedon="inner")
u1 = GridFunction(fes1, "u1")
u1.Set(x*y)
Draw(u1)

[6]:

BaseWebGuiScene

Note how \(u_1\) is displayed only in the inner region.

Do not be confused with ndof for spaces on subdomains:

[7]:

fes = H1(mesh)
fes1.ndof, fes.ndof

[7]:

(132, 132)

Although ndofs are the same for spaces on one subdomain and the whole domain, FreeDofs show that the real number of degrees of freedom for the subdomain space is much smaller:

[8]:

print(fes1.FreeDofs())

0: 11110000111110000000000000000000000000000000000001
50: 10000000000000000000000000000000000000000000000000
100: 00000000000000000000000000000000

One can also glean this information by examining CouplingType of each degrees of freedom of the space, which reveals that there are many unknowns of type COUPLING_TYPE.UNUSED_DOF.

[9]:

for i in range(fes1.ndof):
    print(fes1.CouplingType(i))

COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.WIREBASKET_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF
COUPLING_TYPE.UNUSED_DOF

It is possible to remove these dofs (thus making GridFunction vectors smaller) as follows.

[10]:

fescomp = Compress(fes1)
print(fescomp.ndof)

Integrating on regions¶

You have already seen how boundary regions are used in setting Dirichlet boundary conditions.

[11]:

fes = H1(mesh, order=3, dirichlet="b|l|r")
u, v = fes.TnT()
gfu = GridFunction(fes)

Such boundary regions or subdomains can also serve as domains of integration.

[12]:

term1 = u1 * v * dx(definedon=mesh.Materials("inner"))
term2 = 0.1 * v * ds(definedon=mesh.Boundaries("t"))
f = LinearForm(term1 + term2).Assemble()

Here the functional \(f\) is defined as a sum of two integrals, one over the inner subdomain, and another over the top boundary:

\[f(v) = \int_{\Omega_{inner}} u_1 v \; dx + \int_{\Gamma_{top}} \frac{v}{10}\; ds\]

The python objects dx and ds represent volume and boundary integration, respectively. They admit a keyword argument definedon that is either a Region or a str, as can be seen from the documentation: look for definedon: Optional[Union[ngsolve.comp.Region, str]] in the help output below.

[13]:

help(dx)

Help on DifferentialSymbol in module ngsolve.comp object:

class DifferentialSymbol(pybind11_builtins.pybind11_object)
 |  Method resolution order:
 |      DifferentialSymbol
 |      pybind11_builtins.pybind11_object
 |      builtins.object
 |
 |  Methods defined here:
 |
 |  __call__(...)
 |      __call__(self: ngsolve.comp.DifferentialSymbol, definedon: Optional[Union[ngsolve.comp.Region, str]] = None, element_boundary: bool = False, element_vb: ngsolve.comp.VorB = <VorB.VOL: 0>, skeleton: bool = False, bonus_intorder: int = 0, intrules: dict[ngsolve.fem.ET, ngsolve.fem.IntegrationRule] = {}, deformation: ngsolve.comp.GridFunction = None, definedonelements: pyngcore.pyngcore.BitArray = None) -> ngsolve.comp.DifferentialSymbol
 |
 |  __init__(...)
 |      __init__(self: ngsolve.comp.DifferentialSymbol, arg0: ngsolve.comp.VorB) -> None
 |
 |  ----------------------------------------------------------------------
 |  Static methods inherited from pybind11_builtins.pybind11_object:
 |
 |  __new__(*args, **kwargs) class method of pybind11_builtins.pybind11_object
 |      Create and return a new object.  See help(type) for accurate signature.

Thus - dx(definedon=mesh.Materials("inner")) may be replaced by dx('inner'),

and similarly,

ds(definedon=mesh.Boundaries("t") may be replaced by ds('t').

[14]:

f = LinearForm(u1*v*dx('inner') + 0.1*v*ds('t')).Assemble()

[15]:

a = BilinearForm(fes)
a += grad(u)*grad(v)*dx
a.Assemble()
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw(gfu);

Operations with `Region`s¶

Query mesh for all regions:

[16]:

mesh.GetBoundaries()   # list boundary/interface regions

[16]:

('interface', 'interface', 'interface', 'interface', 'b', 'l', 'r', 't')

[17]:

mesh.GetMaterials()    # list all subdomains

[17]:

('inner', 'outer')

Print region information:

[18]:

print(mesh.Materials("inner").Mask())
print(mesh.Materials("[a-z]*").Mask())  # can use regexp
print(mesh.Boundaries('t|l').Mask())

0: 10
0: 11
0: 00000101

Add regions:

[19]:

io = mesh.Materials("inner") + mesh.Materials("outer")
print(io.Mask())

0: 11

Take complement of a region:

[20]:

c = ~mesh.Materials("inner")
print(c.Mask())

0: 01

Subtract regions:

[21]:

diff = mesh.Materials("inner|outer") - mesh.Materials("outer")
print(diff.Mask())

0: 10

Set a piecewise constant CoefficientFunction using the subdomains:

[22]:

domain_values = {'inner': 3.7,  'outer': 1}
cf = mesh.MaterialCF(domain_values)
Draw(cf, mesh);

Use boundary regions to define a coefficient function:

[23]:

# Make a linear function equaling 2 at the bottom right vertex
# of (0,2) x (0,2) and equals zero at the remaining vertices:
bdry_values = {'b': x, 'r': 2-y}
cf = mesh.BoundaryCF(bdry_values, default=0)

Drawing such a BoundaryCF does not show anything useful. However, we can use a GridFunction and Set, which then lets us view an extension of these boundary values into the domain.

[24]:

g = GridFunction(H1(mesh), name='bdry')
g.Set(cf, definedon=mesh.Boundaries('b|r'))
Draw(g);

1.5 Spaces and forms on subdomains¶

Naming subdomains and boundaries¶

A finite element space on a subdomain¶

Integrating on regions¶

Operations with Regions¶

Operations with `Region`s¶