Multiple Domains General Question

Hi PaulR,

I think you can assemble different equations on different domains once you have mesh function on cells of the mesh, i.e., the labels on the cells. Then you can use dx(label) to tell the compiler to integrate only on cells with that label. I do not have complete demo, but i am sure this works.

commented Jan 29, 2015 by jying FEniCS User (2,020 points)

Thank you for your comment. After a lot of trial and error I finally have some working code with two subdomains and two submeshes.

But how do I create and solve one system of equations from both subdomains and how do I set the interface conditions?

Could you please just point me into the right direction again?

Here is my code:

import math as m
import numpy as np
from dolfin import *

################################################################################
# set measurements and create mesh
thickness = 0.1
length = 0.3

Lx = length
Ly = 0.01
Lz = 0.01
Nx = 360
Ny = 1
Nz = 1
mesh = BoxMesh(0., 0., 0., Lx, Ly, Lz, Nx, Ny, Nz)

# simulation name
simname = 'test'

################################################################################
# set periodic boundary conditions and create constrained function space
class PeriodicDomain(SubDomain):
    # set periodic boundary conditions

    def inside(self, x, on_boundary):
        # return True if on left or bottom boundary AND NOT on one of the two
        # slave edges
        return bool((near(x[1], 0) or near(x[2], 0.)) and
                    (not ((near(x[1], Ly) and near(x[2], 0.)) or
                     (near(x[1], 0) and near(x[2], Lz)))) and on_boundary)

    def map(self, x, y):
        if near(x[1], Ly) and near(x[2], Lz):
            y[0] = x[0]
            y[1] = x[1] - Ly
            y[2] = x[2] - Lz
        elif near(x[1], Ly):
            y[0] = x[0]
            y[1] = x[1] - Ly
            y[2] = x[2]
        elif near(x[2], Lz):
            y[0] = x[0]
            y[1] = x[1]
            y[2] = x[2] - Lz
        else:
            y[0] = -1000
            y[1] = -1000
            y[2] = -1000

periodicdomain = PeriodicDomain()

# ################################################################################

class Air(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= thickness

class Absorber(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] <= thickness

class ABC_Sommerfeld(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[0], Lx)

air = Air()
absorber = Absorber()
abcs = ABC_Sommerfeld()

# create sumdomains
domains = CellFunction('size_t', mesh)
domains.set_all(0)
absorber.mark(domains, 1)
air.mark(domains, 2)

# create submeshes
abs_mesh = SubMesh(mesh, domains, 1)
air_mesh = SubMesh(mesh, domains, 2)

# create subdomains of submeshes
abssubdomains = CellFunction('size_t', abs_mesh)
abssubdomains.set_all(0)
dx_abssubdomain = Measure('dx')[abssubdomains]

airsubdomains = CellFunction('size_t', air_mesh)
airsubdomains.set_all(0)
dx_airsubdomain = Measure('dx')[airsubdomains]

air_boundaries = FacetFunction('size_t', air_mesh)
air_boundaries.set_all(0)
abcs.mark(air_boundaries, 1)
ds_air = Measure('ds')[air_boundaries]

# create two independent function spaces
W_air = FunctionSpace(air_mesh, 'Lagrange', 2, constrained_domain=periodicdomain)
W2_air = W_air * W_air
W_abs = FunctionSpace(abs_mesh, 'Lagrange', 2, constrained_domain=periodicdomain)
W2_abs = W_abs * W_abs

################################################################################
# Define constants and expressions for variational problem
c = 343.4
rho = 1.204
xi = 15000.
kappa = 1.5
sigma = 0.99

zr = Constant(kappa)
zi = Constant(-xi*sigma/(rho))

################################################################################
frequencies = np.array([3000.])
# frequencies = np.arange(20, 5624, 20)
for frequency in frequencies:
    print frequency

    omega = 2.*m.pi*frequency
    agemo = Constant(1./omega)
    k = omega/c

    p0r = 1.
    p0i = 1.
    f1 = Expression('2*k_*p0i_', k_=k, p0i_=p0i)
    f2 = Expression('2*k_*p0r_', k_=k, p0r_=p0r)
    ksquared = Constant(k**2.)
    k = Constant(k)

    (pr_abs, pi_abs) = TrialFunction(W2_abs)
    (pr_air, pi_air) = TrialFunction(W2_air)
    (w1_abs, w2_abs) = TestFunction(W2_abs)
    (w1_air, w2_air) = TestFunction(W2_air)

    # absorber
    a1 = (pr_abs*zr*w1_abs*ksquared - zi*agemo*pi_abs*w1_abs*ksquared + 
         zi*agemo*pr_abs*w2_abs*ksquared + zr*pi_abs*w2_abs*ksquared -
         inner(nabla_grad(pr_abs), nabla_grad(w1_abs)) - 
         inner(nabla_grad(pi_abs), nabla_grad(w2_abs)))*dx_abssubdomain(0)

    # air
    a2 = (pr_air*w1_air*ksquared  + pi_air*w2_air*ksquared -
         inner(nabla_grad(pr_air), nabla_grad(w1_air)) - 
         inner(nabla_grad(pi_air), nabla_grad(w2_air)))*dx_airsubdomain(0)

    a3 = (k*pi_air*w1_air - k*pr_air*w2_air)*ds_air(1) 
    a = a2 + a3
    L1 = f1*w1_air*ds_air(1) - f2*w2_air*ds_air(1)

    A_air = assemble(a, exterior_facet_domains=air_boundaries)
#     A_abs = assemble(a1)

    b_air = assemble(L1, exterior_facet_domains=air_boundaries)

    # Compute solution
    p_air = Function(W2_air)
    P_air = p_air.vector()
    solve(A_air, P_air, b_air)
    pr_air, pi_air = p_air.split()

#     p_abs = Function(W2_abs)
#     P_abs = p_abs.vector()
#     solve(A_abs, P_abs)
#     pr_abs, pi_abs = p_abs.split()

    plot(pr_air)
    plot(pi_air)
#     plot(pr_abs)
#     plot(pi_abs)
    interactive()

commented Jan 30, 2015 by PaulR FEniCS Novice (420 points)

Do you mean you have more than one equations on each subdomains? If this is the case, I am not exactly sure (did not test this before), but you can use MixedFunctionSpace to declare the function space you need and then proceed as the normal case by using each function component corresponding to different equations. Interface conditions for each equations then can be coupled into equations just during the formulations of variational problem. I am also learning using FEniCS. Sorry if this cannot help you.

commented Jan 30, 2015 by jying FEniCS User (2,020 points)