Hi

I recently installed the mixed-dimensional branch and I am trying to make sense out of it. I want to solve a problem on a simple 2D geometry. The geometry is a rectangle fixed at the left and right edges under distributed pressure load at the top edge. This rectangle is divided equally into 2 rectangles (left and right subdomains). The left subdomain is a linear elastic material but the right subdomain is a hyperelastic material . This figure shows the problem that I am looking at:

In other words, I want to solve a linear and a nonlinear problem simultaneously. I have already defined the geometry mesh, marked domains and boundaries and also defined the variational problem on each subdomain but in the last step, I am not sure how to solve it (Maybe because there are new syntax that I am not familiar with). Does anybody know how I can solve it?

Here is the minimal work:

```
from dolfin import *
mesh = RectangleMesh(Point(0., 0.), Point(10, 2), 20, 5)
n = FacetNormal(mesh)
# Ditributed load at the top
p = Constant(10)
# Define boundaries
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 0.0)
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[0], 10.0)
class Top_Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 2.0) and x[0] >= 0.0 and x[0] <= 5.0
class Top_Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[1], 2.0) and x[0] >= 5.0 and x[0] <= 10.0
boundaries = MeshFunction('size_t', mesh, mesh.topology().dim()-1)
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim())
left = Left()
right = Right()
top_left = Top_Left()
top_right = Top_Right()
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)
top_left.mark(boundaries, 3)
top_right.mark(boundaries, 4)
# Define subdomains
class LEFT_HALF(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 0.0 and x[0] <= 5.0 and x[1] >= 0.0 and x[1] <= 2.0 else False
class RIGHT_HALF(SubDomain):
def inside(self, x, on_boundary):
return True if x[0] >= 5.0 and x[0] <= 10.0 and x[1] >= 0.0 and x[1] <= 2.0 else False
left_half = LEFT_HALF()
right_half = RIGHT_HALF()
# Marking left and right subdomains
subdomains.set_all(0)
left_half.mark(subdomains, 1)
right_half.mark(subdomains, 2)
# Linear elasticity on the left domain
def eps(v):
return sym(grad(v))
E = Constant(1e5)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
# Function spaces for the left and right subdomains
V_left = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
V_right = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
W = MixedFunctionSpace(V_left, V_right)
# test and trial functions for the left and right subdomains
du_left, du_right = TrialFunctions(W)
v_left, v_right = TestFunctions(W)
# BCs
bc_left = DirichletBC(W.sub_space(0), Constant((0.,0.)), boundaries , 1)
bc_right = DirichletBC(W.sub_space(1), Constant((0.,0.)), boundaries , 2)
bcs = [bc_left,bc_right]
# Variational form for the left subdomain
a = inner(sigma(du_left), eps(v_left))*dx(1) + inner(sigma(du_left), eps(v_left))*dx(1) # dx(1) corresponds to the left subdomain
l = inner(-p*n, v_left) * ds(3) # ds(3) corresponds to the top edge of the left subdomain
u_mixed = Function(W)
#solve(a == l, u_mixed, bcs)
# Hyperelasticity on the right subdomain
I = Identity(2) # Identity tensor
F = I + grad(u_mixed) # Deformation gradient
C = F.T*F # Right Cauchy-Green tensor
# Invariants of deformation tensors
Ic = tr(C)
J = det(F)
# Stored strain energy density (compressible neo-Hookean model)
psi = (mu/2)*(Ic - 3) - mu*ln(J) + (lmbda/2)*(ln(J))**2
# Total potential energy ### dx(2) corresponds to the right subdomain
Pi = psi*dx(2) - inner(-p*n, v_right) * ds(4) # ds(4) corresponds to the top edge of the right subdomain
# directional derivative
F = derivative(Pi, u_mixed, v_right)
# Compute Jacobian of F
J = derivative(F, u_mixed, du_right)
# Solve variational problem
#solve(F == 0, u_mixed, bcs, J=J)
# How to solve it?!
```