from dolfin import *
import numpy as np
import fenics as fe
#continue the code
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[0] <= length/2 + DOLFIN_EPS
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[0] > length/2 - DOLFIN_EPS
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
materials = MeshFunction('size_t', mesh, mesh.topology().dim())
materials.set_all(0)
subdomain_0.mark(materials, 1)
subdomain_1.mark(materials, 2)
# Boundary conditions
class Boundaryleft(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 0, DOLFIN_EPS)
class RightBoundary(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and near(x[0], 1, DOLFIN_EPS)
V = VectorFunctionSpace(mesh, 'Lagrange', 2)
# ADDED
facets = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
facets.set_all(0)
Boundaryleft().mark(facets, 1)
RightBoundary().mark(facets, 2)
bc1 = DirichletBC(V.sub(0), Constant(0), facets, 1)
# Define measure for boundary integration
# applying neumann boundary condition on the right side
neumann = Expression("100000*x[0]", degree=1)
ds = Measure('ds', domain=mesh, subdomain_data=facets)
dx = Measure('dx', domain=mesh, subdomain_data=materials)
#applying direchla boundary contion on the left side
bc = bc1
# rest of the code ....
u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, 0))
T = Constant((0, 0, 0))
a = inner(sigma(u, E_global, nu), epsilon(v)) * dx
n = FacetNormal(mesh) # Normal vector to the boundary
L = dot(f, v) * dx + neumann*dot(n, v)*ds(2)
# Compute solution
u_solution = Function(V)
solve(a == L, u_solution, bc)
this is what i did, one facet is neumann and the other is direchlet