I am using FEniCs to apply a load on particular point of a square mesh. When checking people has proposed to use Pointwise but it seems to be a method to fix a displacement, others has proposed to use PointSource but this also is used for one component of a vector or I don’t how to use it correctly.
Any help please.
Here is the code I am using:
from future import print_function
#!/usr/bin/env python3
-- coding: utf-8 --
“”"
Created on Thu Dec 24 18:36:09 2020
@author: pierre
“”"
#""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""#
from fenics import *
from mshr import *
import numpy as np
Nx = 1
Ny = 1
L = 1Nx
H = 1Ny
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny)
E = Constant(1e9)
nu = Constant(0.3)
beta = E/(1-nu**2)
C = beta*as_matrix([[1, nu, 0],[nu, 1, 0],[0, 0, (1-nu)/2]])
#""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""#
def eps(v):
return sym(grad(v))
def strain2voigt(e):
“”“e is a 2nd-order tensor, returns its Voigt vectorial representation”""
return as_vector([e[0,0],e[1,1],2*e[0,1]])
def voigt2stress(s):
“”"
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
“”"
return as_tensor([[s[0],s[2]],[s[2],s[1]]])
def sigma(v):
return voigt2stress(dot(C,strain2voigt(eps(v))))
#""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""#
class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],L) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],0) and on_boundary
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],0) and on_boundary
class point(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],Nx) and near(x[1],Ny)
#""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""""#
exterior facets MeshFunction
facets = MeshFunction(“size_t”, mesh, 1)
facets.set_all(0)
Top().mark(facets, 1)
Left().mark(facets, 2)
Bottom().mark(facets, 3)
point().mark(facets,4)
ds = Measure(‘ds’)[facets]
#""" “”"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""#
Define function space
V = VectorFunctionSpace(mesh, ‘Lagrange’, 2)
#""" “”"""""""""""""""""""""""""""""""""""""""""""""""""#
Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name=‘Displacement’)
a = inner(sigma(du), eps(u_))*dx
#""" “”""""""""""""""""""""""""""""""""""""""""""""""#
uniform traction on top boundary
T = Constant((10,10))
l = dot(T, u_)*ds(4)
#""" “”""""""""""""""""""""""""""""""""""""""""""""""""""""""#
symmetry boundary conditions
bc = [DirichletBC(V, Constant((0.,0.)), facets, 2),
DirichletBC(V, Constant((0.,0.)), facets, 3),
DirichletBC(V, T, point(), method=“pointwise”)]
#""" “”"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""#
solve(a == l, u, bc)
#"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
import matplotlib.pyplot as plt
plt.figure(1)
plot(u, title = “u”)
plt.show()
