Fix material property in partial domain

1

Hello,

I am working on an optimization problem, which finds the optimized material property while part of the material property is fixed. See following for an example code with poisson’s equation and some illustrations.

Screen Shot 2022-09-20 at 10.55.45 AM
Screen Shot 2022-09-20 at 10.55.45 AM833×702 105 KB

In the above drawing, kappa2 is meant to be an area with a constant value, while kappa1 takes the form of any function input. Note that using an expression with a conditional statement would not work properly in my optimization implementation.
One way I’m picturing the implementation be something like this: kappa = kappa1+kappa2, where kappa1 only defines in the purple subdomain, and kappa2 only defines in the yellow subdomain.
I am struggling how to make this implementation work, any advice is appreciated.

from dolfin import *
from mshr import*

domain = Rectangle(Point(0., 0.), Point(1., 1.))
domain.set_subdomain(1, Rectangle(Point(0.6, 0.1), Point(0.8, 0.3)) )

mesh = generate_mesh(domain, 30)
mf = MeshFunction("size_t", mesh, 2, mesh.domains())
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(6.4*2,4.8*2),dpi=800)
plot(mesh, "2D mesh")
plot(mf, "Subdomains")
plt.show()


V = FunctionSpace(mesh, "Lagrange", 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS


u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)


u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)

kappa = Function(V)

a = inner( kappa*  grad(u), grad(v))*dx
L = f*v*dx + g*v*ds

Forward = a-L
2

You could use sub-domains to mark each region, and only integrate with a constant kappa over one, and a kappa defined as a function on the other.

3

Hi Jorgen, thank you for your respond, do you have any example I could follow? I am clueless on how to start this.

4

Consider:

from IPython import embed
import matplotlib.pyplot as plt
from mshr import*
from dolfin import *

domain = Rectangle(Point(0., 0.), Point(1., 1.))
domain.set_subdomain(1, Rectangle(Point(0.6, 0.1), Point(0.8, 0.3)))

mesh = generate_mesh(domain, 30)
mf = MeshFunction("size_t", mesh, 2, mesh.domains())

V = FunctionSpace(mesh, "Lagrange", 1)


def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS


u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)


u = TrialFunction(V)
v = TestFunction(V)
f = Expression(
    "10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
g = Expression("sin(5*x[0])", degree=2)

kappa0 = Function(V)
kappa1 = Constant(0.1)
dx = Measure("dx", domain=mesh, subdomain_data=mf)
a = inner(kappa0 * grad(u), grad(v))*dx(0) + \
    inner(kappa1 * grad(u), grad(v))*dx(1)
L = f*v*dx + g*v*ds

Forward = a-L
print(assemble(1*dx(0)), assemble(1*dx(1)))

where I’ve printed the areas of each of the integration domains at the end

0.9600000000000003 0.04000000000000004
5

Is it possible to define kappa0 only in dx(0)?

6

It might be possible by using MeshView. However, if you want to use tools for automatic differentiation such as dolfin_adjoint the method described above should give the correct sensitivities ( sensitivities of kappa0 in dx(1) will be Zero).

7

I will give that a try. Also when I tried to plot the solution u, I received this error:

Traceback (most recent call last):
  File "demo_poisson.py", line 57, in <module>
    plot(u)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 438, in plot
    return _plot_matplotlib(object, mesh, kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 304, in _plot_matplotlib
    return mplot_function(ax, obj, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/dolfin/common/plotting.py", line 152, in mplot_function
    return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 321, in tricontourf
    return TriContourSet(ax, *args, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 40, in __init__
    ContourSet.__init__(self, ax, *args, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/contour.py", line 816, in __init__
    kwargs = self._process_args(*args, **kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 52, in _process_args
    tri, z = self._contour_args(args, kwargs)
  File "/opt/anaconda3/envs/fenicsproject/lib/python3.8/site-packages/matplotlib/tri/tricontour.py", line 99, in _contour_args
    raise ValueError('z array must not contain non-finite values '
ValueError: z array must not contain non-finite values within the triangulation

my plotting code is simply:

import matplotlib.pyplot as plt
plt.figure()
plot(u)
plt.show()

Is there a specific way to plot u when it is solved with multiple subdomain?

8

Please provide a minimal example that reproduces your error

9

sorry I found a bug in my testing code, plotting problem solved now, I will try on passing this code to differentiation as you suggested, thank you!