Dear @dokken
Thanks a lot for the reply and guidance. I have built a minimal working code that describes the issue for which I have requested help. The code is attached below:
I have a function u defined on the function space, Vsub in the top submesh. I have obtained this function by solving a PDE in this submesh.
from dolfin import *
from fenics import *
from ufl import nabla_div
from ufl import sinh
r = Constant(25.2)
g = -100
C = Constant(19.34)
kappa = 0.04
# create mesh
mesh = RectangleMesh(Point(0,0), Point(100e-6, 200e-6), 100,100)
regions = MeshFunction('size_t', mesh, mesh.topology().dim())
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim() - 1)
class A(SubDomain):
def inside(self, x, on_boundary):
return x[1] < 100e-6 #Bottom half of domain
a = A()
class B(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 100e-6 #Top half of domain
b = B()
class TopBoundary(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 200e-6)
top = TopBoundary()
class MiddleBoundary(SubDomain):
def inside(self, x, on_boundary):
return near (x[1], 100e-6)
middle = MiddleBoundary()
regions.set_all(0)
b.mark(regions, 1) #Top region is marked 1; bottom remains marked as 0
boundaries.set_all(0)
middle.mark(boundaries, 2) #Middle boundary in the domain
top.mark(boundaries, 1) #Top most boudnary in the domain
# Working on the top submesh - To solve a PDE
# define a submesh, composed of region 1 - top region
new_region = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region.set_all(0)
new_region.array()[regions.array() == 1] = 1
submesh = SubMesh(mesh, new_region, 1)
# define new meshfunctions on this submesh
submesh_regions = MeshFunction('size_t', submesh, submesh.topology().dim())
submesh_regions.set_all(0)
submesh_boundaries = MeshFunction('size_t', submesh, submesh.topology().dim() - 1)
submesh_boundaries.set_all(0)
middle.mark(submesh_boundaries, 2)
top.mark(submesh_boundaries, 1)
ds = Measure("ds")[submesh_boundaries]
Vsub = FunctionSpace(submesh, 'CG', 1)
u = TrialFunction(Vsub)
v = TestFunction(Vsub)
u = Function(Vsub)
dx = Measure("dx")[submesh_regions]
bcs = []
F = kappa*inner(nabla_grad(u), nabla_grad(v))*dx + g*v*ds(1) + r*sinh(C*u)*v*ds(2)
J = derivative(F, u)
problem = NonlinearVariationalProblem(F, u, bcs, J)
solver = NonlinearVariationalSolver(problem)
solver.solve()
#Plotting solution u obtained on Vsub - defined on top submesh
import matplotlib.pyplot as plt
from matplotlib import cm
p = plot(u, cmap = cm.inferno)
plt.colorbar(p)
plt.show()
# Working on the bottom submesh
# define a submesh, composed of region 0 - bottom region
new_region2 = MeshFunction('size_t', mesh, mesh.topology().dim())
new_region2.set_all(0)
new_region2.array()[regions.array() == 0] = 1
submesh2 = SubMesh(mesh, new_region2, 1)
# define new meshfunctions on this submesh
submesh_regions2 = MeshFunction('size_t', submesh2, submesh2.topology().dim())
submesh_regions2.set_all(0)
submesh_boundaries2 = MeshFunction('size_t', submesh2, submesh2.topology().dim() - 1)
submesh_boundaries2.set_all(0)
middle.mark(submesh_boundaries2, 2)
ds2 = Measure("ds")[submesh_boundaries2]
Vsub2 = FunctionSpace(submesh2, 'CG', 1)
u2 = TrialFunction(Vsub2)
v2 = TestFunction(Vsub2)
dx2 = Measure("dx")[submesh_regions2]
u2 = Function(Vsub2)
Following this, I need to define a function u2 in the bottom submesh that takes the values of u at the shared boundary and zero everywhere else. Something like this:
values of u2(at the shared boundary) = values of u at the shared boundary (This u is defined in the top submesh)
Since, both the boundaries are marked and share the same co-ordinates, I thought of extracting the dofs and modifying the function u2 as follows:
from numpy import where
top_dofs = where(submesh_boundaries.array()==2)
bottom_dofs = where(submesh_boundaries2.array()==2)
u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]
It shows the following error.
u2.vector()[bottom_dofs[0]] = u.vector()[top_dofs[0]]
IndexError: Vector index out of range
I’m doing some fundamental mistake in assigning the function values based on the dofs, and not able to update the function u2 properly.
Please help me resolving this issue. Thanks in advance !