# 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 0x7fd5800635f0>


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


## 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.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
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)

11


## 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) from pybind11_builtins.pybind11_type
|      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.Assemble()
gfu.vec.data = a.mat.Inverse(fes.FreeDofs()) * f.vec
Draw(gfu);


## Operations with Regions¶

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

0: 10
0: 11
0: 00000101


[19]:

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

0: 11


Take complement of a region:

[20]:

c = ~mesh.Materials("inner")

0: 01


Subtract regions:

[21]:

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

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);