Hi. I am Solving the elasticity problem in a 2d model to represents the in situ stress in a underground mining. In the model , the right and left boundaries have a different maximum height as you can see in the picture. I want to apply a compressure condition in this axis as a function depend from x[1] but in my unexperience y havent found the correc way to do this.
My script is:
from dolfin import *
mesh = Mesh("test_anisotropic5.xml")
#plot(mesh)
cd=MeshFunction('size_t',mesh,"test_anisotropic5_physical_region.xml")
#plot(cd)
#fd=MeshFunction('size_t',mesh,"test_anisotropic5_facet_region.xml")
#plot(fd)
class Young(UserExpression): # UserExpression instead of Expression
def __init__(self, markers, **kwargs):
super().__init__(**kwargs) # This part is new!
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = 30e3
elif self.markers[cell.index] ==2:
values[0] = 15e3
elif self.markers[cell.index] ==3:
values[0] = 15e3
elif self.markers[cell.index] ==4:
values[0] = 30e3
elif self.markers[cell.index] ==5:
values[0] = 15e3
elif self.markers[cell.index] ==6:
values[0] = 15e3
else:
values[0] = 30e3
#Inputs Values
E=Young(cd,degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)
x_max = Constant(914.5695) #End x axis point
yl_max= Constant(562) #End y axis point in x=0
yr_max= Constant(672.322) #End y axis point in x=x_max
k = Constant(1.5)
rho_g = 0.027
f = Constant((0, -rho_g))
left_function = Expression("yl_max-0.027*k*x[1]",0) #Compression in left boundary
rigth_function = Expression("yr_max-0.027*k*x[1]",0) #Compression in rigth Boundary
rigth = AutoSubDomain(lambda x: near(x[0],0.0))
left = AutoSubDomain(lambda x: near(x[0],x_max))
# Definition of Neumann boundary condition domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
rigth.mark(boundaries, 2)
ds = ds(subdomain_data=boundaries)
def bottom(x, on_boundary):
return near(x[1], 0.)
bc = DirichletBC(V, Constant((0.,0.)), bottom)
def eps(v):
return sym(grad(v))
def sigma(v):
return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx+(inner(left_function,u_))*ds(1) #i need to aplicate in the rigth axis too
l = inner(f, u_)*dx
This is what im trying to implementate the neumann BC, but im not sure that is the correct wat. I hope you can help me…
You should define the compression on the left and right boundaries as vector-valued Expressions, e.g.
left_function = Expression(("0.027*k*(yl_max-x[1])", 0), degree=1) #Compression in left boundary
rigth_function = Expression(("-0.027*k*(yr_max-x[1])", 0), degree=1) #Compression in rigth Boundary
(Note that I have amended the expressions to produce something like the pressure distributions shown in your schematic.) Since the expressions you have provided are linear, you should specify at least a degree >= 1 space for the expressions to be interpolated into.
I would refer to the boundary near x_0=0 as the left boundary, and the boundary near x_0=x_{max} as the right boundary, but that’s just my personal choice. Regardless, you may want to double-check your subdomain definitions to ensure they are consistent with the compressions you have defined.
Finally, the Neumann boundary should be incorporated in the linear form l, not in the bilinear form a:
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx + inner(left_function,u_)*ds(1) + inner(rigth_function, u_)*ds(2)
2 Likes
Hi Conpierce, thanks for your answer, i understood this perfect!..
But i still have a problem when i run my scrip whit te funcions i see this error:
INFO:dijitso:------------------- Start compiler output ------------------------
/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp: In member function ‘virtual void dolfin::dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:
/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp:61:29: error: ‘k’ was not declared in this scope
61 | values[0] = 0.027*k*(yl_max-x[1]);
| ^
/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp:61:32: error: ‘yl_max’ was not declared in this scope
61 | values[0] = 0.027*k*(yl_max-x[1]);
| ^~~~~~
INFO:dijitso:/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp: In member function ‘virtual void dolfin::dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2::eval(Eigen::Ref<Eigen::Matrix<double, -1, 1> >, Eigen::Ref<const Eigen::Matrix<double, -1, 1> >) const’:
/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp:61:29: error: ‘k’ was not declared in this scope
61 | values[0] = 0.027*k*(yl_max-x[1]);
| ^
/tmp/tmprozac7kq/dolfin_expression_7f5e1531f48f6d4818a1106f1cae47d2.cpp:61:32: error: ‘yl_max’ was not declared in this scope
61 | values[0] = 0.027*k*(yl_max-x[1]);
| ^~~~~~
I im tryng to wrhite the constant values k and yl_max whit “Constant”, but it still dont work… I hope you can help me…
def bottom(x, on_boundary):
return near(x[1], 0.)
bc1 = DirichletBC(V, Constant((0.,0.)), bottom)
bc=[bc1]
x_max = 914.5695 #End x axis point
k = Constant(1.5)
yl_max= 562.0 #End y axis point in x=0
yr_max= 672.322
left_function = Expression(("0.027*k*(yl_max-x[1])", 0), degree=1) #Compression in left boundary
rigth_function = Expression(("-0.027*k*(yr_max-x[1])", 0), degree=1) #Compression in rigth Boundary
rigth = AutoSubDomain(lambda x: near(x[0],0.0))
left = AutoSubDomain(lambda x: near(x[0],x_max))
# Definition of Neumann boundary condition domain
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
left.mark(boundaries, 1)
rigth.mark(boundaries, 2)
ds = ds(subdomain_data=boundaries)
rho_g = 0.027
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx + inner(left_function,u_)*ds(1) + inner(rigth_function, u_)*ds(2)
My mistake. You need to pass yl_max and yr_max as kwargs to Expression, i.e.:
left_function = Expression(("0.027*k*(yl_max-x[1])", 0), degree=1, k=k, yl_max=yl_max)
rigth_function = Expression(("-0.027*k*(yr_max-x[1])", 0), degree=1, k=k, yr_max=yr_max)
2 Likes
Thanks so much now it works perfect!!
1 Like
Hi conpierce. Now i am tryng to implementate the same Neuman bc in a 3D model, which is like a cube whit topografy in the top face like:
In 2D model i define the compress funcion whit a y_max value, but in 3d i have diferents max Z value along the face. I think in some script that can read the max value of Z along the face and define the function depending on these values. I hope you can help me…
