elasticity cantilever beam - constant load at the corner

Asked by Niraj Kumar Jha

Hi teammates,
I have solved a cantilever beam problem fixed at the left end and distributed load at the top... here is the code
from dolfin import *
from dolfin_utils import meshconvert

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create classes for defining parts of the boundaries and the interior
# of the domain

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], 2.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):
        return ((x[1], (1.0)) and (x[0], (2.0)))
        #return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

class Point(SubDomain):
    def inside(self, x, on_boundary):
        return pt

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
pt = Point()

# Create mesh and define function space
#mesh = dolfin.UnitSquare(2, 2)
mesh = Rectangle(0.0, 0.0, 2.0, 1.0, 20, 10)
plot(mesh)
print mesh.cells()
print mesh.data()
myArray= mesh.coordinates()
pt = [myArray[5][0],myArray[5][1]]
print pt
mesh.geometry()

# Initialize mesh function for interior domains
domains = CellFunction("uint", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

# Define function space and basis functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
U = TrialFunction(V)
v = TestFunction(V)

# Define boundary conditions
c = Expression(("0.0", "0.0"))
bc = DirichletBC(V, c, boundaries, 1)

# Define new measures associated with the interior domains and exterior boundaries

dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]

# Define functions

v = TestFunction(V) # Test function
U = TrialFunction(V)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

def epsilon(v):return 0.5*(grad(v)+transpose(grad(v)))
def sigma(v):return 2*mu*epsilon(v)+lmbda*tr(epsilon(v))*Identity(V.cell().d)

# Define variational problem
g_T= Constant((0.0, -0.5))
a = inner(grad(v),sigma(U))*dx(1) # nodes
L= dot(v,g_T)*ds(2)

# Solve problem
U = Function(V)
solve(a == L, U, bc)

# Plot solution and mesh

file = File("elasticity1.pvd")
file << U

# Hold plot
plot(U, interactive=True)

I need to solve the same prb with constant point load acting (in negative y direction only ) @ co-ordinate (2,1)..please suggest me to solve the problem...

thnks

Question information

Language:
English Edit question
Status:
Answered
For:
FEniCS Project Edit question
Assignee:
No assignee Edit question
Last query:
Last reply:
Revision history for this message
Nick Davies (ntd14) said :
#1

You can do it the same way as you implement the force on the top. just define a new boundary at the point you want.
eg
class point(SubDomain):
    def inside(self, x, on_boundary):
        return near(2.0, 1.0)

point.mark(boundaries, 5)

L= dot(v,g_T)*ds(5)

Revision history for this message
Niraj Kumar Jha (niraj-jha) said :
#2

thanks for your reply.
I am able to run this code but there is no result (all the displacements are zero)..can you please give a reason??

from dolfin import *
from dolfin_utils import meshconvert

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
ffc_options = {"optimize": True, \
               "eliminate_zeros": True, \
               "precompute_basis_const": True, \
               "precompute_ip_const": True}

# Create classes for defining parts of the boundaries and the interior
# of the domain

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], 2.0)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 0.0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], 1.0)

class Obstacle(SubDomain):
    def inside(self, x, on_boundary):

        # return on_boundary and (near(x[0], 2.0) or near(x[1], 0.0))
        #return True if x[1] >= 1.0 and x[0] >= 2.0 else False
        return ((x[1], (1.0)) and (x[0], (2.0)))
        #return (between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0)))

class Point(SubDomain):
    def inside(self, x, on_boundary):
        return near(2.0, 1.0)

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
obstacle = Obstacle()
pt = Point()

# Create mesh and define function space
#mesh = dolfin.UnitSquare(2, 2)
mesh = Rectangle(0.0, 0.0, 2.0, 1.0, 20, 10)
plot(mesh)
print mesh.cells()
print mesh.data()
#myArray= mesh.coordinates()
#pt = [myArray[5][0],myArray[5][1]]
#print pt
mesh.geometry()

# Initialize mesh function for interior domains
domains = CellFunction("uint", mesh)
domains.set_all(0)
obstacle.mark(domains, 1)

# Initialize mesh function for boundary domains
boundaries = FacetFunction("uint", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
pt.mark(boundaries, 5)

# Define function space and basis functions
V = VectorFunctionSpace(mesh, "Lagrange", 1)
U = TrialFunction(V)
v = TestFunction(V)

# Define boundary conditions
c = Expression(("0.0", "0.0"))
bc = DirichletBC(V, c, boundaries, 1)

# Define new measures associated with the interior domains and exterior boundaries

dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]

# Define functions

v = TestFunction(V) # Test function
U = TrialFunction(V)

# Elasticity parameters
E, nu = 10.0, 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

def epsilon(v):return 0.5*(grad(v)+transpose(grad(v)))
def sigma(v):return 2*mu*epsilon(v)+lmbda*tr(epsilon(v))*Identity(V.cell().d)

# Define variational problem
g_T= Constant((0.0, -0.5))
a = inner(grad(v),sigma(U))*dx(1) # nodes
L= dot(v,g_T)*ds(5)

# Solve problem
U = Function(V)
solve(a == L, U, bc)

# Plot solution and mesh

file = File("elasticity.pvd")
file << U

# Hold plot
plot(U, interactive=True)

Revision history for this message
Johannes Ring (johannr) said :
#3

FEniCS no longer uses Launchpad for Questions & Answers. Please consult the documentation on the FEniCS web page for where and how to (re)post your question: http://fenicsproject.org/support/

Can you help with this problem?

Provide an answer of your own, or ask Niraj Kumar Jha for more information if necessary.

To post a message you must log in.