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):
return sym(grad(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?
please help me
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,
Based on your advice to see the link below:
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):
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 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.
Please I need your help
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):
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 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:
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
msh = meshio.read("/home/pierre/Desktop/GMSH1/carre.msh")
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)],
cell_data={"name_to_read":[triangle_data]})
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):
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],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
Cell data: name_to_read, 1
I tried to correct it based on the error message but unsuccessfully.
Please need your help again.
Thanks
You need to load the mesh from XDMF. i.e.
from dolfin import *
# First write mesh to xdmf
# ...
# Then define functions and classes
# ...
# Then load mesh
mesh = Mesh()
with XDMFFile("test.xdmf") as infile:
infile.read(mesh)
# 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>.
Please, how to solve this?
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
# ...
# Then load mesh
mesh = Mesh()
with XDMFFile("test.xdmf") as infile:
infile.read(mesh)
# 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