plot sigm.n on the boundary


418
views
1
12 months ago by
HI, I have solved elasticity problem for crack. Now i want to plot σ.n on the boundary But I am getting error: UFLException: Can only integrate scalar expressions. The integrand is a tensor expression with value rank 1 and free indices ().

Kindly find the attached minimal working code:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
set_log_active(False)
#==================================
# mesh (concident nodes)
#===================================
nx = 1
ny = 1
num = 1
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh, 2, 2)
editor.init_vertices(10)
editor.init_cells(8)
editor.add_vertex(0, 0, 0.0)
editor.add_vertex(1, 0.5, 0.0)
editor.add_vertex(2, 0.5, 0.5)
editor.add_vertex(3, 0.0, 0.5)
editor.add_vertex(4, -0.5, 0.5)
editor.add_vertex(5, -0.5, 0.0)
editor.add_vertex(6, -0.5, -0.5)
editor.add_vertex(7, 0.0, -0.5)
editor.add_vertex(8, 0.5, -0.5)
editor.add_vertex(9, -0.5, 0.0)
editor.add_cell(0, 0, 1, 3)
editor.add_cell(1, 1, 2, 3)
editor.add_cell(2, 0, 3, 4)
editor.add_cell(3, 0, 4, 5)
editor.add_cell(4, 0, 9, 7)
editor.add_cell(5, 6, 7, 9)
editor.add_cell(6, 0, 7, 8)
editor.add_cell(7, 0, 8, 1)
editor.close()
num_refine = 3
h=0.5**num_refine
for i in range(num_refine):
mesh=refine(mesh)
print mesh.hmin()
plot(mesh)
print mesh.num_cells()
interactive()

V = VectorFunctionSpace(mesh,'P',1)
E = 10e5
nu = 0.28
tol = 1e-10
mu = E/(2*(1+nu))
lamb = E*nu/((1+nu)*(1-2*nu))

def clamped_boundary(x,on_boundary):
    return abs(x[1]+0.5) < tol and on_boundary

bc1 = DirichletBC(V,Constant((0,0)),clamped_boundary)

def right(x,on_boundary):
    return abs(x[1]-0.5) < tol and on_boundary

bc2 = DirichletBC(V.sub(1),Constant((0.1)),right)
bc=[bc1,bc2]
def epsilon(u):
    return 0.5*(grad(u)+grad(u).T)
def sigma(u):
    return lamb*(tr(grad(u))*Identity(d))+2*mu*epsilon(u)

u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
form = inner(sigma(u),epsilon(v))*dx
a = lhs(form)
L = rhs(form)
u = Function(V)
solve(a == L,u,bc)
u1,u2=u.split()
plot(u,mode='displacement',interactive=True)
n = FacetNormal(mesh)
KK = TensorFunctionSpace(mesh, "Lagrange", 2)
sigma_w = project(sigma(u), KK)
plot(sigma_w[0, 0], key = 'sigma11', title = 'sigma11')
plot(sigma_w[1, 1], key = 'sigma22', title = 'sigma22')
Traction = dot(sigma_w,n)*ds
print "Traction", shape(Traction)
Traction = project(Traction)
plot(Traction)​
Community: FEniCS Project
Can you format the code properly?  Thanks.
written 12 months ago by pf4d  
Thanks for the suggestion. Now, I have formatted it properly.
written 12 months ago by hirshikesh  
Not quite, that for loop needs an indent...
written 12 months ago by pf4d  
What exactly do you want to plot, the traction \( \mathbf{t} = \sigma \mathbf{n} \) on the boundary or its integral? Plotting the integral -- which is what you're trying to do -- doesn't make much sense, since it is just an integral of a vector quantity giving you a constant vector.
written 12 months ago by Adam Janecka  
Multiplication by ds doesn't restrict the quantities to the boundary but integrates over it.
written 12 months ago by Adam Janecka  

2 Answers


0
12 months ago by
pf4d  
If you want a better response to your questions, I suggest stating them as minimal working examples.
For this question, you are asking why the following code fails :

from fenics import *

mesh     = UnitSquareMesh(10,10, 'crossed')
V        = TensorFunctionSpace(mesh, "Lagrange", 2)
n        = FacetNormal(mesh)
tensor   = interpolate(Constant([[0,1],[2,3]]), V)
traction = dot(tensor,n)*ds

This is because you can only integrate scalar expressions, and the integrand of "traction" is a vector.

Instead, you can do something like this :

traction_x = dot(tensor[0,:],n)*ds​

​for each component of your traction.
Thank you for the reply and suggestions, I want to plot σ.n also. When I try to plot I got error: AttributeError: 'Form' object has no attribute 'shape'
written 12 months ago by hirshikesh  
What set of procedures let you to this error?  Code here is helpful.

Here's one way to plot what you want :

from fenics import *

mesh     = UnitSquareMesh(10,10, 'crossed')
V        = TensorFunctionSpace(mesh, "Lagrange", 2)
Q        = FunctionSpace(mesh, "CG", 1)
n        = FacetNormal(mesh)
u        = TestFunction(Q)
tensor   = interpolate(Constant([[0,1],[2,3]]), V)
traction = dot(tensor[0,:],n)*u*ds

# a function to hold the computation :
v        = Function(Q)

# set the local degrees of freedom to the assembled normal component of 
# tensor :
v.vector().set_local(assemble(traction).array())

# plot it :
plot(v, interactive=True)​
written 12 months ago by pf4d  
Hi, Thaks for the reply.

When I add use this:
from dolfin import *
import numpy as np
import matplotlib.pyplot as plt
set_log_active(False)
#==================================
# mesh (concident nodes)
#===================================
nx = 1
ny = 1
num = 1
mesh = Mesh()
editor = MeshEditor()
editor.open(mesh, 2, 2)
editor.init_vertices(10)
editor.init_cells(8)
editor.add_vertex(0, 0, 0.0)
editor.add_vertex(1, 0.5, 0.0)
editor.add_vertex(2, 0.5, 0.5)
editor.add_vertex(3, 0.0, 0.5)
editor.add_vertex(4, -0.5, 0.5)
editor.add_vertex(5, -0.5, 0.0)
editor.add_vertex(6, -0.5, -0.5)
editor.add_vertex(7, 0.0, -0.5)
editor.add_vertex(8, 0.5, -0.5)
editor.add_vertex(9, -0.5, 0.0)
editor.add_cell(0, 0, 1, 3)
editor.add_cell(1, 1, 2, 3)
editor.add_cell(2, 0, 3, 4)
editor.add_cell(3, 0, 4, 5)
editor.add_cell(4, 0, 9, 7)
editor.add_cell(5, 6, 7, 9)
editor.add_cell(6, 0, 7, 8)
editor.add_cell(7, 0, 8, 1)
editor.close()
num_refine = 3
h=0.5**num_refine
for i in range(num_refine):
    mesh=refine(mesh)
print mesh.hmin()
plot(mesh)
print mesh.num_cells()
interactive()

V = VectorFunctionSpace(mesh,'P',1)
E = 10e5
nu = 0.28
tol = 1e-10
mu = E/(2*(1+nu))
lamb = E*nu/((1+nu)*(1-2*nu))

def clamped_boundary(x,on_boundary):
    return abs(x[1]+0.5) < tol and on_boundary

bc1 = DirichletBC(V,Constant((0,0)),clamped_boundary)

def right(x,on_boundary):
    return abs(x[1]-0.5) < tol and on_boundary

bc2 = DirichletBC(V.sub(1),Constant((0.1)),right)
bc=[bc1,bc2]
def epsilon(u):
    return 0.5*(grad(u)+grad(u).T)
def sigma(u):
    return lamb*(tr(grad(u))*Identity(d))+2*mu*epsilon(u)

u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
form = inner(sigma(u),epsilon(v))*dx 
a = lhs(form)
L = rhs(form)
u = Function(V)
solve(a == L,u,bc)
u1,u2=u.split()
plot(u,mode='displacement')
n = FacetNormal(mesh)
KK = TensorFunctionSpace(mesh, "Lagrange", 2)
sigma_w = project(sigma(u), KK)
plot(sigma_w[0, 0], key = 'sigma11', title = 'sigma11')
plot(sigma_w[1, 1], key = 'sigma22', title = 'sigma22')
Traction_x = dot(sigma_w[0,:],n)*ds
Traction_y = dot(sigma_w[1,:],n)*ds
Z = FunctionSpace(mesh, 'CG',1)
zz = Function(Z)

zz.vector().set_local(assemble(Traction_x))

# plot it :
plot(zz, interactive=True)

​

 I have got error: TypeError: (2) numpy array of 'double' expected. Make sure that the numpy array use dtype=float_.

written 12 months ago by hirshikesh  
man what a bummer...
written 12 months ago by pf4d  
-1
3 months ago by
zayyf  
editor.open(mesh, 2, 2)
NotImplementedError: Wrong number or type of arguments for overloaded function 'MeshEditor_open'.
what is the reason?
Please login to add an answer/comment or follow this question.

Similar posts:
Search »