Help a starter!

variational formulation
21

For instance Here I have delta for point force, based on my understanding value is the force x0 is the position of force applied. I am expecting a change in the displacement when I change x0 but there is no change. Is my understanding is correct? if not how can I fix this code?

from dolfin import *
from ufl import nabla_div
import numpy as np
from dolfin import *
# Load mesh and define function space
L =10.
W = 1.
# Mesh for cantilever beam
Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W),100, 3, 3)
#Omega.init()
V = VectorFunctionSpace(Omega, "CG", 1)
tol = 1.e-8 # tolerance
#Dirichlet on left boundary
def l_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)
class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        print("init function")
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs) 
    def eval(self, values, x):
       
        eps = self.eps
        values[0] = 0
        values[1] = 0
        values[2] = -1

    def value_shape(self): 
        return (3, ) 
         

delta = Delta(eps=1E-4, x0=[-0,-0.,-10], degree=3)

zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, l_boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0., -1.0))
# Elasticity parameters
E = 1e5
nu = 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Strain
def eps(u):
    return sym(grad(u))
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(u):
    return lmbda*tr(eps(u))*Identity(3) + 2.0*mu*eps(u)



# Define variational problem

a = inner(sigma(u), eps(v))*dx

l=  dot(f, v)*dx +dot(v,delta)*ds
# Compute solution
u = Function(V)
solve(a == l, u, bc)
print("Maximal deflection:", -u(L,0.5,0.5)[2])
22

A Dirac delta function is a smoothed point load. If you do as I said on January 31st you will get a point load as well.

I the example below, you are not using a delta function

You are setting the function to constant -1 everywhere. Please read through your code carefully.

23

thank you, it would be great if you can show how to do, I could not find how to fix it

24

I am gonna stop replying to you, as you are clearly not putting in any effort to solve this yourself. I have already shown you how to do this in Applying a force on a face

25

Yes, there is this command there for value, for 2d case :
eps/pi/(np.linalg.norm(x-self.x0)2 + eps2)

from dolfin import *
from ufl import nabla_div
import numpy as np
from dolfin import *
# Load mesh and define function space
L =10.
W = 1.
# Mesh for cantilever beam
Omega = BoxMesh(Point(0, 0, 0), Point(L, W, W),100, 3, 3)
#Omega.init()
V = VectorFunctionSpace(Omega, "CG", 1)
tol = 1.e-8 # tolerance
#Dirichlet on left boundary
def l_boundary(x, on_boundary):
    return on_boundary and near(x[0], 0.0)
class Delta(UserExpression):
    def __init__(self, eps, x0, **kwargs):
        print("init function")
        self.eps = eps
        self.x0 = x0
        UserExpression.__init__(self, **kwargs) 
    def eval(self, values, x):
       
        eps = self.eps
        values[0] =-0
        #eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
        values[1] = 0
        values[2] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)

    def value_shape(self): 
        print("ah")
        return (3, ) 
         

delta = Delta(eps=1E-4, x0=[-10,-0.,-0], degree=3)

zero = Constant((0.0, 0.0, 0.0))
bc = DirichletBC(V, zero, l_boundary)
u = TrialFunction(V)
v = TestFunction(V)
f = Constant((0.0, 0., -1.0))
# Elasticity parameters
E = 1e5
nu = 0.3
mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))
# Strain
def eps(u):
    return sym(grad(u))
model = "plane_stress"
mu = E/2/(1+nu)
lmbda = E*nu/(1+nu)/(1-2*nu)
if model == "plane_stress":
    lmbda = 2*mu*lmbda/(lmbda+2*mu)
def sigma(u):
    return lmbda*tr(eps(u))*Identity(3) + 2.0*mu*eps(u)



# Define variational problem

a = inner(sigma(u), eps(v))*dx

l=  dot(f, v)*dx +dot(v,delta)*ds
# Compute solution
u = Function(V)
solve(a == l, u, bc)
print("Maximal deflection:", -u(L,0.5,0.5)[2])

but still displacement not affected much by the point (x0 )

26

Increase the magnitude of the dirac delta function then. A point source in 3D is less effective than 2D

27

if I want to apply a specific value for the force is there a way to do it? what is the the magnitude of force here