how can i create a big domain and an obstacle inside and solve a problem only inside of the obstacle?

dear team,

i have a question :( i'm new to fenics and dolfin-adjoint.

the following code is an example i found from page: http://www.dolfin-adjoint.org/en/latest/documentation/optimisation.html

and from the following page i found an example of multiple domains: http://fenicsproject.org/documentation/dolfin/1.2.0/python/demo/pde/subdomains-poisson/python/documentation.html

i want to combine the problems. it means, i want create a great domain und solve the problem of the heat equation inside an obstacle of this great domain. anywhere must be a mistake in the code and my thinking. because it must be the same result. but it isnt. and also my u_desired is "nan" (not a number). what is the mistake? and is it possible to draw a line or is it possible to make a distinction of the two domains in the plots? please help me to fix my problem! please scroll down and there is the code i've written.

best regards

christian student

#######################################################################################
# Start of Code
#######################################################################################

#######################################################################################
# import all what is possible and of interest maybe
from fenics import *
from dolfin import *
from dolfin_adjoint import *
from math import *
# another imports
import sys
import scipy
import numpy
#######################################################################################

#######################################################################################
# subdomain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], -3.0)

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 3.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], -3.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 3.0)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return (between(x[0], (-1.0, 1.0)) and between(x[1], (-1.0, 1.0)))
#######################################################################################

#######################################################################################
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
#######################################################################################

#######################################################################################
# define mesh
n = 600
mesh = RectangleMesh(Point(-3.0, -3.0), Point(3.0, 3.0), n, n)
#mesh = UnitSquareMesh(600, 600)
#######################################################################################

#######################################################################################
# Initialize mesh function for interior domains
domains = CellFunction("size_t", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)
#######################################################################################

#######################################################################################
# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
#######################################################################################

#######################################################################################
# Define input data
a0 = Constant(1.0)
a1 = Constant(1.0)
f = Constant(0)
#######################################################################################

#######################################################################################
# define functions and function spaces
V = FunctionSpace(mesh, "CG", 1)
W = FunctionSpace(mesh, "DG", 1)
u = Function(V, name = "State")
m = Function(W, name = "Control")
v = TestFunction(V)
#######################################################################################

#######################################################################################
# dirichlet boundary conditions
#bc = DirichletBC(V, 0.0, "on_boundary")
bc = [DirichletBC(V, 0.0, boundaries, 1),
      DirichletBC(V, 0.0, boundaries, 2),
      DirichletBC(V, 0.0, boundaries, 3),
      DirichletBC(V, 0.0, boundaries, 4)]
#######################################################################################

#######################################################################################
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure("dx")[domains]
#######################################################################################

#######################################################################################
# Run the forward model once to create the simulation record
F = (inner(a0*grad(u), grad(v)) - m*v)*dx(1) + (inner(a1*grad(u), grad(v)) - f*v)*dx(0)
solve(F == 0, u, bc)
#######################################################################################

#######################################################################################
# The functional of interest is the normed difference between desired
# and simulated temperature profile
x = triangle.x
u_desired = Expression("exp(-1/(1-x[0]*x[0])-1/(1-x[1]*x[1]))")
J = Functional((0.5*inner(u-u_desired, u-u_desired))*dx(1)*dt[FINISH_TIME])
#######################################################################################

#######################################################################################
# Run the optimisation
J_Tilde = ReducedFunctional(J, Control(m, value = m))
# Make sure you have scipy >= 0.11 installed
m_opt = minimize(J_Tilde, method = "L-BFGS-B", tol = 2e-12, bounds = (0.0, 0.5), 
                 options = {"disp": True})
# replace 
F_opt = replace(F, {m: m_opt})
#a, L = lhs(F_opt), rhs(F_opt)
solve(F_opt == 0, u, bc)
#solve(a == L, u, bc)
#######################################################################################

#######################################################################################
# plots
plot(m_opt, title = "control")
plot(mesh)
u_d_plot = interpolate(u_desired, V)
plot(u_d_plot, title = "u_desired")
plot(u, title = "u_real")
# option for plots
interactive()
#######################################################################################

#######################################################################################
# end of code
#######################################################################################

Comments (4)