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

`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

``````E.vector()[:] = 2
E.vector()[pos]=5
``````

would assign the value 5 the cell with index `pos` and 2 to all other cells

Hi Dokken,
Thank you very much.

Hi,
I am getting an error in a readable file that I am not understanding.
Code

``````import meshio
for cell in msh.cells:
triangle_cells = cell.data
for key in msh.cell_data_dict["gmsh:physical"].keys():
if key == "triangle":
triangle_data = msh.cell_data_dict["gmsh:physical"][key]
mesh =meshio.Mesh(points=msh.points,
cells=[("triangle", triangle_cells)],
meshio.write("test.xdmf",mesh)

E = 2.1e11    #p[0]
nu = 0.3    #p[1]
beta = (E)/(1-(nu*nu))
C = as_matrix([[beta, beta*nu, 0],[beta*nu, beta, 0],[0, 0, beta*(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],10) 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
#********************************************************************#
# exterior facets MeshFunction
facets = MeshFunction("size_t", mesh, 1)
``````

Error message

``````
File "/home/pierre/Desktop/MesEssais/essai28_jun/essai28_jun22.py", line 59, in <module>
facets = MeshFunction("size_t", mesh, 1)

File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/mesh/meshfunction.py", line 18, in __new__
return fn(mesh, dim)

TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
1. dolfin.cpp.mesh.MeshFunctionSizet(arg0: dolfin.cpp.mesh.Mesh, arg1: int)
2. dolfin.cpp.mesh.MeshFunctionSizet(arg0: dolfin.cpp.mesh.Mesh, arg1: int, arg2: int)
3. dolfin.cpp.mesh.MeshFunctionSizet(arg0: dolfin.cpp.mesh.Mesh, arg1: str)
4. dolfin.cpp.mesh.MeshFunctionSizet(arg0: dolfin.cpp.mesh.Mesh, arg1: int, arg2: dolfin.cpp.mesh.MeshDomains)
5. dolfin.cpp.mesh.MeshFunctionSizet(arg0: dolfin.cpp.mesh.Mesh, arg1: dolfin::MeshValueCollection<unsigned long>)

Invoked with: <meshio mesh object>
Number of points: 146
Number of cells:
triangle: 250
``````

I tried to correct it based on the error message but unsuccessfully.

Thanks

You need to load the mesh from `XDMF`. i.e.

``````from dolfin import *
# First write mesh to xdmf
# ...
# Then define functions and classes
# ...

mesh = Mesh()
with XDMFFile("test.xdmf") as infile:
# Then create mesh function
facets = MeshFunction("size_t", mesh, 1)
``````

Thank you dokken,
I have followed your instructions, but I have another error message:

``````Shapes do not match: <ListTensor id=140219301787456> and <Sym id=140219301787200>.
Traceback (most recent call last):

File "/home/pierre/Desktop/MesEssais/essai28_jun/essai28_jun22.py", line 86, in <module>
a = inner(sigma(du), eps(u_))*dx

File "/usr/lib/python3/dist-packages/ufl/operators.py", line 158, in inner
return Inner(a, b)

File "/usr/lib/python3/dist-packages/ufl/tensoralgebra.py", line 147, in __new__
error("Shapes do not match: %s and %s." % (ufl_err_str(a), ufl_err_str(b)))

File "/usr/lib/python3/dist-packages/ufl/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))

UFLException: Shapes do not match: <ListTensor id=140219301787456> and <Sym id=140219301787200>.
``````

You need to provide a minimal code example that reproduces the error, as this is related to the variational form

Ok, based on your previous advice, I have developed the code bellow:

``````from dolfin import *
# First write mesh to xdmf
# ...
# Then define functions and classes
# ...

mesh = Mesh()
with XDMFFile("test.xdmf") as infile:

# 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)
ds = Measure("ds", subdomain_data = facets)

# Define function space
V = VectorFunctionSpace(mesh, 'CG', 1)

# 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
Ttrac = Constant((0., 1e7))
l = dot(Ttrac, u_)*ds(1)
``````

It is the code that gives the previous error I sent you.

If your mesh is 2 dimensional, you Need to prune the z-coordinates of points before writing to file (you can use ` mesh.prune_z_0()` for the mesh above.
You can verify this by printing `mesh.geometry().dim()` after loading the mesh into dolfin; it should be 2

1 Like