# How to apply a load on a single point or a node

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.
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 = 1
Ny

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):
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()

Please format your code appropriately, using 3x` encapsulation.you should also consider the answers in: Implement point source in the variational formulation

Ok, Thank you I am going to consider your instructions.

Hi Dokken,
I considered the answers like you advise me. The fragment of code that I modified to use is the one below. I thought that X0 is the place where I can apply the load but it seems not to be so. Even if I change the values of X0, there is no effect on the point where I want to apply the load.

``````class Delta(UserExpression):
def __init__(self, eps, x0, **kwargs):
self.eps = eps
self.x0 = x0
UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
eps = self.eps
values[0] = 10    #eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values[1] = 10
# values[2] = 0

def value_shape(self): return (2, )

delta = Delta(eps=1E-4, x0=np.array([1, 1]), degree=5)
# l = inner(Constant(100)*delta, u_)*dx
# delta = Constant((10, 10))

l = inner(delta, u_)*dx
``````

Just have a look at the plot of a source function, in for instance paraview:

``````from dolfin import *
import numpy as np
mesh = UnitSquareMesh(20, 20)
V = VectorFunctionSpace(mesh, "CG", 1)

class Delta(UserExpression):
def __init__(self, eps, x0, **kwargs):
self.eps = eps
self.x0 = x0
UserExpression.__init__(self, **kwargs)
def eval(self, values, x):
eps = self.eps
values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values[1] = 0

def value_shape(self): return (2, )

delta = Delta(eps=1E-4, x0=np.array([1, 1]), degree=5)

u = project(delta, V)
File("u.pvd") << u
``````

As you are not specifying what is not working I can’t you help you more.

1 Like

Thanks for your spontaneous answer. To be more clear, my idea is to apply a load in x and y direction at the point of coordinates (1, 1). So I have changed the lines below

``````values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values[1] = 0
``````

into

``````values[0] = 10    #eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
values[1] = 10
``````

but as the last two lines are independent of x0, I am not getting a satisfied result.
Or may be you can clearly explain me how the line below is working in the code

``````values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
``````

Hope that I have made sens to my question now.
Thanks

As you observe in the code I gave you, it gave you a point source in one direction.
The x dependence above is what creates the point source as you divide epsilon by the distance from x0. If you apply values[0]=10, values[1]=10, you set a constant source all over the domain

Hi dokken,
Excuse me for I come back to the same question. I want to control the value of the values[0], so I print it in order to control the value I will obtain but I obtained too many values.

`````` values[0] = eps/pi/(np.linalg.norm(x-self.x0)**2 + eps**2)
print(values[0])
``````

See the result I obtain

``````
3.183098830006919e-05
3183.0988618379065
4.9735918939093545e-05
8.841941037273602e-05
0.00019894366643089014
0.000795774516515848
2.4867959663827704e-05
4.420970580039168e-05
9.947183632393973e-05
0.000397887307993825
1.940913928310224e-05
2.340513851788801e-05
2.7440507193081033e-05
3.060671953106911e-05
3.183098830006919e-05
3.9788735275614644e-05
4.681027669158277e-05
6.121343847354747e-05
7.957746955651093e-05
0.00015915493513414867
1.591549422961206e-05
3.183098830006919e-05
3183.0988618379065
4.9735918939093545e-05
8.841941037273602e-05
0.00019894366643089014
0.000795774516515848
2.4867959663827704e-05
4.420970580039168e-05
9.947183632393973e-05
0.000397887307993825
1.940913928310224e-05
2.340513851788801e-05
2.7440507193081033e-05
3.060671953106911e-05
3.183098830006919e-05
3.9788735275614644e-05
4.681027669158277e-05
6.121343847354747e-05
7.957746955651093e-05
0.00015915493513414867
1.591549422961206e-05
3.183098830006919e-05
3183.0988618379065
4.9735918939093545e-05
8.841941037273602e-05
0.00019894366643089014
0.000795774516515848
2.4867959663827704e-05
4.420970580039168e-05
9.947183632393973e-05
0.000397887307993825
1.940913928310224e-05
2.340513851788801e-05
2.7440507193081033e-05
3.060671953106911e-05
3.183098830006919e-05
3.9788735275614644e-05
4.681027669158277e-05
6.121343847354747e-05
7.957746955651093e-05
0.00015915493513414867
1.591549422961206e-05
3.183098830006919e-05
3183.0988618379065
4.9735918939093545e-05
8.841941037273602e-05
0.00019894366643089014
0.000795774516515848
2.4867959663827704e-05
4.420970580039168e-05
9.947183632393973e-05
0.000397887307993825
1.940913928310224e-05
2.340513851788801e-05
2.7440507193081033e-05
3.060671953106911e-05
3.183098830006919e-05
3.9788735275614644e-05
4.681027669158277e-05
6.121343847354747e-05
7.957746955651093e-05
0.00015915493513414867
``````

I want to know exactly the value I am giving to values[0]. Is it possible ?
That is I want it to take only one value the I can change as I wish.
Thanks.

`values[0]` is of the dimension of a scalar function space (print x to see this). This means that you are evaluating your function at a set of coordinates at once. This, if you want to control a specific value in this array, you need to identify which entry you would like to change (by inspecting x).

Hi,
Hope you are ok.

I am trying to solve a 3D (Orthotropic) linear elasticity in FEniCs but still getting an error message that I am not understanding very well.
Code

``````from ufl import *
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt

#*********************************************************************#
N = 2
mesh = UnitCubeMesh(N, N, N)
# plt.figure(1)
# plot(mesh)
# plt.show()

#************************************************************#
Ex, Ey, Ez, nuxy,nuxz,nuyz, Gxy, Gxz, Gyz = 100., 10.,1., 0.3, 0.2, 0.1, 5., 4., 3.
S = as_matrix([[1./Ex,-nuxy/Ex,-nuxz/Ex,0.,0.,0.],[-nuxy/Ex,1./Ey,-nuyz/Ey,0.,0.,0.],[-nuxz/Ex,-nuyz/Ey,1/Ez,0.,0.,0.],[0.,0.,0.,1/Gxy,0.,0.],[0.,0.,0.,0.,1/Gxz,0.],[0.,0.,0.,0.,0.,1/Gyz]])
C = inv(S)
print('S',S)
print('C',C)

#******************************************************************#
def eps(v):
def strain3voigt(e):
"""e is a 2nd-order tensor, returns its Voigt vectorial representation"""
return as_vector([e[0,0],e[1,1],e[2,2],2*e[0,1],2*e[0,2],2*e[1,2]])
def voigt3stress(s):
"""
s is a stress-like vector (no 2 factor on last component)
returns its tensorial representation
"""
return as_tensor([[s[0],s[3],0.],[s[3],s[1],0.],[0.,0.,s[2]]])
def sigma(v):
return voigt3stress(dot(C,strain3voigt(eps(v))))

#***************************************************************************#

class Top(SubDomain):
def inside(self, x, on_boundary):
return near(x[2],1) and on_boundary
class Left(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],0) and on_boundary
class Bottom(SubDomain):
def inside(self, x, on_boundary):
return near(x[2],0) and on_boundary
class Right(SubDomain):
def inside(self, x, on_boundary):
return near(x[1],1) and on_boundary
class Direct(SubDomain):
def inside(self, x, on_boundary):
return near(x[0],0) and on_boundary
#***********************************************************************#
# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 2)
facets.set_all(0)
Top().mark(facets, 1)
Left().mark(facets, 2)
Bottom().mark(facets, 3)
Right().mark(facets,4)
Direct().mark(facets,5)
ds = Measure('ds', subdomain_data=facets)
#***********************************************************************#
# Define function space
V = VectorFunctionSpace(mesh, 'CG', 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((0.,0.,1e-3))
l = dot(T, u_)*ds(1)
#*********************************************************************#
# symmetry boundary conditions
bcs = [DirichletBC(V.sub(0), Constant(0.), facets, 2),
DirichletBC(V.sub(1), Constant(0.), facets, 3),
DirichletBC(V.sub(2), Constant(0.), facets, 5)]
#********************************************************************#

u = Function(V)
problem = LinearVariationalProblem(a, l, u, bcs)
solver = LinearVariationalSolver(problem)
solver.solve()
solver.parameters['linear_solver'] = 'cg'
solver.parameters['preconditioner'] = 'ilu'
cg_prm = solver.parameters['krylov_solver']
cg_prm['absolute_tolerance'] = 1E-7
cg_prm['relative_tolerance'] = 1E-4
cg_prm['maximum_iterations'] = 1000
``````

Error message

``````UFLException: adj_expr not implemented for dimension 6.
``````

Is it my Voigt notation that is not correct ?
How to correct my program to solve a 3D orthotropic linear elasticity in FEniCs?

Hi, you cannot use the UFL `inv` operator to invert the 6x6 compliance matrix. As said in the corresponding doc, `inv` is implemented only for 2x2 or 3x3 since it uses a symbolic formula for the inversion. You must therefore either implement yourself the expression of C or go through numpy for instance to invert S as a numpy array first.

1 Like

Hi Bleyerj, Thank you

Hi,
I hope you are ok,
I am trying to solve a linear elasticity problem.

But my objective is to assign different Young modulus to the different parts of the mesh above. For instance I want to assign E1 to the yellow part, E2 to the green part…
Please can you show me the way to do so?
I am reading documents but still not understanding very well.
Thanks

See for instance:

Hi Dokken,
https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomain
I have developed the code below:

``````from __future__ import print_function

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
"""
Created on Thu Jun  3 08:36:47 2021

@author: pierre
"""

#***********************************************************#
from fenics import *
from mshr import *
import numpy as np

#***********************************************************#
# Create the mesh
Nx, Ny = 1 , 1
L, H = 1, 1
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, 'left')
#************************************************************#

E0 = Constant(2.1e11)
nu0 = Constant(0.3)
beta0 = E0/(1-(nu0*nu0))
C0 = as_matrix([[beta0, beta0*nu0, 0],[beta0*nu0, beta0, 0],[0, 0, beta0*(1-nu0)/2]])
#*************************************************************#
E1 = Constant(70e9)
nu1 = Constant(0.3)
beta1 = E1/(1-(nu1*nu1))
C1 = as_matrix([[beta1, beta1*nu1, 0],[beta1*nu1, beta1, 0],[0, 0, beta1*(1-nu1)/2]])

#*************************************************************#

def eps(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 sigma0(v):
return voigt2stress(dot(C0,strain2voigt(eps(v))))

#*************************************************************#

def sigma1(v):
return voigt2stress(dot(C1,strain2voigt(eps(v))))
#*************************************************************#

# Create classes for defining parts of the boundaries
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], 1.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)

#*************************************************************#
# Create classes for defining parts of the interior of the domain
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol

#*************************************************************#

# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
Omega_0 = Omega_0()
Omega_1 = Omega_1()

#*************************************************************#

# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Omega_0.mark(domains, 1)
Omega_0.mark(domains, 2)

#*************************************************************#
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)

#*************************************************************#
# Define function space
V = VectorFunctionSpace(mesh, 'CG', 1)
#*************************************************************#
# Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name = 'Displacement')
#*************************************************************#
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)

#*************************************************************#
# Define Dirichlet boundary conditions at bottom boundary
bcs = DirichletBC(V, Constant((0., 0.)), boundaries, 4)

#*************************************************************#
a = inner(sigma0(du), eps(u_))*dx(1) + inner(sigma1(du), eps(u_))*dx(2)

#*************************************************************#
# uniform traction on top boundary
T = Constant((0, 1e7))
l = dot(T, u_)*ds(2)
#*************************************************************#
solve(a == l, u, bcs)
``````

But I have the following error message:

``````*** Error:   Unable to set given (local) rows to identity matrix.
*** Reason:  some diagonal elements not preallocated (try assembler option keep_diagonal).
*** Where:   This error was encountered inside PETScMatrix.cpp.
*** Process: 0
``````

I tried to solve the problem but unsuccessfully.

Sorry Dokken,
I see my error
Thanks

Hi,
Thank you again for the help,
I have now solved the problem of defining subdomains with different Young modulus.

Now suppose I have a domain meshed with 3 node triangles like the figure below:

And I want to give two different Young modulus to triangles among the triangles in the domain. How can I do that?
I have tried with what is done in the link below:
https://fenicsproject.discourse.group/t/how-to-assign-different-material-properties-to-elements-linear-elasticity/3742
But I was not able to solve my problem.
Hope I make sens.

You need to make a minimal example, as I have already linked you to how to make a discontinuous dolfin function that can store different material properties on each cell. Please follow the structure of the posts highlighted in Read before posting: How do I get my question answered? - #4 by nate

Sorry Dokken,
As you have linked me to how to overcome my problem, I have developed the code below

``````from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
#***********************************************************#
# Create the mesh
Nx, Ny = 4 , 4
L, H = 1, 1
mesh = RectangleMesh(Point(0., 0.), Point(L, H), Nx, Ny, 'left')
#************************************************************#
E0 = Constant(2.1e11)
nu0 = Constant(0.3)
beta0 = E0/(1-(nu0*nu0))
C0 = as_matrix([[beta0, beta0*nu0, 0],[beta0*nu0, beta0, 0],[0, 0, beta0*(1-nu0)/2]])
#*************************************************************#
E1 = Constant(2.1e9)
nu1 = Constant(0.3)
beta1 = E1/(1-(nu1*nu1))
C1 = as_matrix([[beta1, beta1*nu1, 0],[beta1*nu1, beta1, 0],[0, 0, beta1*(1-nu1)/2]])
#*************************************************************#
def eps(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 sigma0(v):
return voigt2stress(dot(C0,strain2voigt(eps(v))))
#*************************************************************#
def sigma1(v):
return voigt2stress(dot(C1,strain2voigt(eps(v))))
#*************************************************************#
# Create classes for defining parts of the boundaries
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], 1.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)
#*************************************************************#
# Create classes for defining parts of the interior of the domain
tol = 1E-14
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol
#*************************************************************#
# Initialize sub-domain instances
left = Left()
top = Top()
right = Right()
bottom = Bottom()
Omega_0 = Omega_0()
Omega_1 = Omega_1()
#*************************************************************#
# Initialize mesh function for interior domains
domains = MeshFunction("size_t", mesh, mesh.topology().dim())
domains.set_all(0)
Omega_1.mark(domains, 1)
#*************************************************************#
# Initialize mesh function for boundary domains
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim()-1)
boundaries.set_all(0)
left.mark(boundaries, 1)
top.mark(boundaries, 2)
right.mark(boundaries, 3)
bottom.mark(boundaries, 4)
#*************************************************************#
# Define function space
V = VectorFunctionSpace(mesh, 'CG', 1)
#*************************************************************#
# Define variational problem
du = TrialFunction(V)
u_ = TestFunction(V)
u = Function(V, name = 'Displacement')
#*************************************************************#
# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure('dx', domain=mesh, subdomain_data=domains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
#*************************************************************#
# Define Dirichlet boundary conditions at bottom boundary
bcs = DirichletBC(V, Constant((0., 0.)), boundaries, 4)
#*************************************************************#
a = inner(sigma0(du), eps(u_))*dx(0) + inner(sigma1(du), eps(u_))*dx(1)
#*************************************************************#
# uniform traction on top boundary
T = Constant((0, -1e7))
l = dot(T, u_)*ds(2)
#*************************************************************#
solve(a == l, u, bcs)
``````

and it is working successfully based on what I have seen here:

https://fenicsproject.org/olddocs/dolfin/latest/python/demos/subdomains-poisson/documentation.html?highlight=subdomain

But my problem now is to assign a different Young modulus to a 3 node triangular element as it is done in the following code:

``````from dolfin import *
E_VALUES = [7.5E6, 12E6, 83E6, 42E6,33E6,98E6]
mesh = BoxMesh(Point(0, 0, 0), Point(1, 1, 1), 1, 1, 1)

DG = FunctionSpace(mesh, "DG", 0)
E = Function(DG)
E.vector()[:] = E_VALUES
``````

The code above is assigning different young modulus to all the elements in the mesh.

Now suppose I have 1000 elements and I want to assign different Young modulus to one of them (E1 to one of them and E2 to the 999 other elements); How can I modify the code above to solve my problem?

Thanks again