Save strees in pvd format

1

I am trying to save my model result in pvd format. The model is about the elasticity and y need plot the displacement a stress. When i save the displacements like:
#u = Function(V, name=“Displacement”)
#solve(a == l, u, bc)
#vtkfile=File(‘solution_displa.pvd’)
#vtkfile << u
Its work. But, y continue the proces to get the principal stress values like:

#stress = sigma(u)
#sx = -stress[0,0]
#sy = -stress[1,1]
#txy = -stress[0,1]
#def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)2 + tau_xy2))

#def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)2 + tau_xy2));
#S1 = sigma1(sx,sy,txy);
#S3 = sigma3(sx,sy,txy)

And
#vtkfile=File("S1.pvd’)
#vtkfile << S1 Dont work.
How i can do this?? Thanks.

The complete code is:

from dolfin import *
mesh = Mesh(“p1_2vetas.xml”)
#plot(mesh)

cd=MeshFunction(‘size_t’,mesh,“p1_2vetas_physical_region.xml”)
plot(cd)
def eps(v):
return sym(grad(v))

class Young(UserExpression): # UserExpression instead of Expression
def init(self, markers, **kwargs):
super().init(**kwargs) # This part is new!
self.markers = markers
def eval_cell(self, values, x, cell):
if self.markers[cell.index] == 1:
values[0] = 30e9
elif self.markers[cell.index] ==2:
values[0] = 15e9
elif self.markers[cell.index] ==3:
values[0] = 15e9
elif self.markers[cell.index] ==4:
values[0] = 30e9
elif self.markers[cell.index] ==5:
values[0] = 15e9
elif self.markers[cell.index] ==6:
values[0] = 15e9
else:
values[0] = 30e9
E=Young(cd,degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =Enu/(1+nu)/(1-2nu)

def sigma(v):
return lmbdatr(eps(v))Identity(2) + 2.0mueps(v)

rho_g = 1e-3
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, ‘Lagrange’, degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx;

def bottom(x, on_boundary):
return near(x[1], 0.)
def left(x, on_boundary):
return near(x[0], 0.)
def right(x, on_boundary):
return near(x[0], 914.5695)
bc1 = DirichletBC(V, Constant((0.,0.)), bottom)
bc2 = DirichletBC(V.sub(0), Constant(0.), left)
bc3 = DirichletBC(V.sub(0), Constant(0.), right)
bc=[bc1,bc2,bc3]

u = Function(V, name=“Displacement”)
solve(a == l, u, bc)
#vtkfile=File(‘solution_displa.pvd’)
#vtkfile << u
p = plot(u,mode=“displacement”,cmap=“jet”)
plt.colorbar(p);
stress = sigma(u)
sx = -stress[0,0]
sy = -stress[1,1]
txy = -stress[0,1]
def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)2 + tau_xy2))

def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)2 + tau_xy2));

S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy);

p = plot(S1,cmap=“jet”);
plt.colorbar(p,label=“σ1(MPa)”);

2

See: Why use DG space to project stress/strain？ - #3 by bhaveshshrimali
or: How to store diferential operator in a Function object - #2 by dokken
Which both suggest projecting the stresses into a suitable function space.

Please use markdown formatting when posting code; i.e.

```python
def f(x):
    return x
```
8

Hi Dokken, thanks for your help and sorry for my beginner error. I alredy get how save the stress values from my model to see it in paraview, however, i dont know how save the principal stress vector from the stress tensor or eigenvalues. I can plot this in an analytical way to visualize in fenics, but i need to do this in paraview, please help me.

from dolfin import *
mesh = Mesh("p1_2vetas.xml")
#plot(mesh)

cd=MeshFunction('size_t',mesh,"p1_2vetas_physical_region.xml")
plot(cd)
def eps(v):
  return sym(grad(v))

class Young(UserExpression): # UserExpression instead of Expression
    def __init__(self, markers, **kwargs):
        super().__init__(**kwargs) # This part is new!
        self.markers = markers
    def eval_cell(self, values, x, cell):
        if self.markers[cell.index] == 1:
            values[0] = 30e3
        elif self.markers[cell.index] ==2:
            values[0] = 15e3
        elif self.markers[cell.index] ==3:
            values[0] = 15e3
        elif self.markers[cell.index] ==4:
            values[0] = 30e3
        elif self.markers[cell.index] ==5:
            values[0] = 15e3
        elif self.markers[cell.index] ==6:
            values[0] = 15e3
        else:
            values[0] = 30e3
E=Young(cd,degree=0)
nu = Constant(0.3)
mu = E/2/(1+nu)
lmbda =E*nu/(1+nu)/(1-2*nu)

def sigma(v):
  return lmbda*tr(eps(v))*Identity(2) + 2.0*mu*eps(v)

rho_g = 0.027
f = Constant((0, -rho_g))
V = VectorFunctionSpace(mesh, 'Lagrange', degree=2)
du = TrialFunction(V)
u_ = TestFunction(V)
a = inner(sigma(du), eps(u_))*dx
l = inner(f, u_)*dx;

def bottom(x, on_boundary):
  return near(x[1], 0.)
def left(x, on_boundary):
  return near(x[0], 0.)
def right(x, on_boundary):
  return near(x[0], 914.5695)    
bc1 = DirichletBC(V, Constant((0.,0.)), bottom)
bc2 = DirichletBC(V.sub(0), Constant(0.), left)
bc3 = DirichletBC(V.sub(0), Constant(0.), right)
bc=[bc1,bc2,bc3]

u = Function(V, name="Displacement")
solve(a == l, u, bc)

Vsig = TensorFunctionSpace(mesh, "DG", degree=0)
sig = Function(Vsig, name="Stress")
sig.assign(project(-sigma(u), Vsig));
sx = sig[0,0]
sy = sig[1,1]
txy = sig[0,1]
def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))
def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2));
S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy)

file_results = XDMFFile("Stress_cartesiano.xdmf")
file_results.parameters["flush_output"] = True
file_results.parameters["functions_share_mesh"] = True
file_results.write(u, 0.)
file_results.write(sig, 0.)
9

Hi @ncarmona,

I am struggling with the same issue. Have you found a solution to this problem?

10

Try this:

#Principal Stress

Vsig = TensorFunctionSpace(mesh, "DG",degree=1)
sig = Function(Vsig, name="Cauchy_Stress")
sig.assign(project(-sigma(u), Vsig)); 

sx = sig.sub(0)
sy = sig.sub(3)
txy = sig.sub(1)

def sigma1(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 + sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2))

def sigma3(sigma_x, sigma_y, tau_xy): # magnitude of first principal stress
    return((sigma_x+sigma_y)/2 - sqrt(((sigma_x-sigma_y)/2)**2 + tau_xy**2));

S1 = sigma1(sx,sy,txy);
S3 = sigma3(sx,sy,txy)

Vp = TensorFunctionSpace(mesh, "DG",degree=1)
sp = Function(Vsig, name="Principal_Stress")
sp.assign(project(as_tensor(((S1,0.),(0.,S3))), Vp));