Appling load inside the domain

General
1

Hello everyone,

I am trying to apply a load on the edge within the domain as shown in the attached figure. Outer boundary of the green rectangle is constrained.

q1
q11160×1160 2.38 KB

My code is as follows

from fenics import *
mesh = Mesh('domain_nm_b.xml')
V = VectorFunctionSpace(mesh, 'CG', 1)
tol = 1E-7
def boundary_D_1(x, on_boundary):
    return x[0] <= 0.0 + tol
def boundary_D_2(x, on_boundary):
    return x[0] >= 4000.0 - tol
def boundary_D_3(x, on_boundary):
    return x[1] <= 0.0 + tol
def boundary_D_4(x, on_boundary):
    return x[1] >= 4000.0 - tol    
bc_1 = DirichletBC(V, Constant((0.,0.)), boundary_D_1,  method='pointwise')
bc_2 = DirichletBC(V, Constant((0.,0.)), boundary_D_2,  method='pointwise')
bc_3 = DirichletBC(V, Constant((0.,0.)), boundary_D_3,  method='pointwise')
bc_4 = DirichletBC(V, Constant((0.,0.)), boundary_D_4,  method='pointwise')
bc_disp = [bc_1, bc_2, bc_3, bc_4]

class boundary_N_y(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 3000.-tol and x[0] <= 3000. + tol
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
B_N_Y = boundary_N_y()
B_N_Y.mark(boundary_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)
E = 210000.
nu = 0.3
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))     
u, du, v = Function(V), TrialFunction(V), TestFunction(V)	
T = Constant((300., 0))
L = dot(T, u)*ds(1)
def eps(u):
    return sym(grad(u))

def sigma(u):
    return 2.0*mu*(eps(u)) + lmbda*tr(eps(u))*Identity(2)
elastic_energy = 0.5*inner(sigma(u), eps(u))*dx - L

E_u = derivative(elastic_energy,u,v)
Jd = derivative(E_u, u, du)

problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters
solver_u.solve()

this code is not applying the load on the required edge and solver is giving the following response:

 Newton iteration 0: r (abs) = 0.000e+00 (tol = 1.000e-10) r (rel) = -nan (tol = 1.000e-09)
  Newton solver finished in 0 iterations and 0 linear solver iterations.

Mesh can be taken from here. Please help me to fix this issue.

Thank you,

Manish

2

You need to replace the “ds” input in Measure to “dS”, as you want to apply the load on an interior facet (a facet connected to two cells).

3

Thank you @dokken . I noticed that along with this, one need to be careful regarding the side on which side of facet force need to apply. I am attaching the updated code to apply force on the right side of the facet for the future readers.

from fenics import *
mesh = Mesh('domain_nm_b.xml')
V = VectorFunctionSpace(mesh, 'CG', 1)
tol = 1E-7
def boundary_D_1(x, on_boundary):
    return x[0] <= 0.0 + tol
def boundary_D_2(x, on_boundary):
    return x[0] >= 4000.0 - tol
def boundary_D_3(x, on_boundary):
    return x[1] <= 0.0 + tol
def boundary_D_4(x, on_boundary):
    return x[1] >= 4000.0 - tol    
bc_1 = DirichletBC(V, Constant((0.,0.)), boundary_D_1,  method='pointwise')
bc_2 = DirichletBC(V, Constant((0.,0.)), boundary_D_2,  method='pointwise')
bc_3 = DirichletBC(V, Constant((0.,0.)), boundary_D_3,  method='pointwise')
bc_4 = DirichletBC(V, Constant((0.,0.)), boundary_D_4,  method='pointwise')
bc_disp = [bc_1, bc_2, bc_3, bc_4]

class boundary_N_y(SubDomain):
    def inside(self, x, on_boundary):
        return x[0] >= 3000.-tol and x[0] <= 3000. + tol
boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
B_N_Y = boundary_N_y()
B_N_Y.mark(boundary_markers, 1)

dS_ = Measure('dS', domain=mesh, subdomain_data=boundary_markers)

E = 210000.
nu = 0.3
lmbda = E * nu / ((1 + nu) * (1 - 2 * nu))
mu = E / (2 * (1 + nu))     
u, du, v = Function(V), TrialFunction(V), TestFunction(V)	
T = Constant((300., 0))
L = dot(T, u('+'))*dS_(1) #dot(T, u('-'))*dS_(1) ##for the left side use u('-')
def eps(u):
    return sym(grad(u))

def sigma(u):
    return 2.0*mu*(eps(u)) + lmbda*tr(eps(u))*Identity(2)
elastic_energy = 0.5*inner(sigma(u), eps(u))*dx - L

E_u = derivative(elastic_energy,u,v)
Jd = derivative(E_u, u, du)

problem_u = NonlinearVariationalProblem(E_u, u, bc_disp, Jd)
solver_u = NonlinearVariationalSolver(problem_u)
prm = solver_u.parameters
solver_u.solve()

Thank you,

Manish