Exporting stress and time to pvd file leads to separate colour palette for each time step?


374
views
0
12 months ago by
Hi, I am interested in stress (von_Mises) and displacement (u) of a solid over time, but if I save them like so 
ufile << (u, t)
mfile << (von_Mises, t)​

I can play an animation of u in Paraview, but it's impossible with von_Mises becauses it seems like every time step receives its own colour palette. I found an intermediary solution by running 

sed -i 's/f_[0-9]*/f_41/g' mises*.vtu

where f_41 is the name of the colour palette in the first vtu file. Doing that I can "play" an animation of the stress over time in Paraview, but the displayed colour just stays the same and the range of values changes instead. Why does that happen when exporting a derived quantity and is there a way to fix it? 

Here's a full example script:

from fenics import *
from dolfin import *
from mshr import *

from matplotlib import pyplot as plt
import numpy as np
from scipy.interpolate import interp1d
rom subprocess import run


# Units
cm = 1e-2
um = 1e-4 * cm
dyn = 1
pa = 10 * dyn/cm**2
s = 1

# Scaled variables
r0 = 20*um
R = r0/r0
Le = 10*R
W = 0.2*R
Disp = 4*um/r0
w_ast = 10*um/r0
ast0 = Le/2-w_ast/2
gap = 1*um/r0
tf = 50*s
asta1 = ast0
asta2 = ast0+w_ast
Y = 1.0e6 * pa
nu = 0.49
lam = Y*nu/((1+nu)*(1-2*nu))
mu = Y/(2*(1+nu))
Mu = mu/mu
Lam = lam/mu

tf = 50
nt = 10
dt = tf/nt
t = 0.0

def disp(t):
    return np.sin(t)*R

# Create mesh
N = 512
deg = 2
elem = "Lagrange"
sim = "test"
geom = Rectangle(Point(0, 0), Point(W, Le))
mesh = generate_mesh(geom, N)
print(mesh)

run("mkdir data/%s" % sim, shell=True)

V = VectorFunctionSpace(mesh, elem, deg)

# Define boundaries
def astr_a(x, on_boundary):
    return near(x[0], W) and (x[1] < asta2 and x[1] > asta1)

def fixed_boundary(x, on_boundary):
    return near(x[1], 0) or near(x[1], Le)

disp_a = Expression(('d', '0'), d=disp(t), degree=deg)
bc1 = DirichletBC(V, Constant((0, 0)), fixed_boundary)
bc2 = DirichletBC(V, disp_a, astr_a, method='pointwise')
bcs = [bc1, bc2]

# Stress and strain
def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return Lam*nabla_div(u)*Identity(d) + 2*Mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
f = Constant((0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx

# Create VTK file for saving solution
ufile = File('data/%s/u_%d.pvd' % (sim, N))
mfile = File('data/%s/mises_%d.pvd' % (sim, N))

# Compute solution
u = Function(V)
t = 0.0
for i in range(nt):
    disp_a.d = disp(i*dt)
    solve(a == L, u, bcs)
    t += dt
    
    # Calculate stress
    W = TensorFunctionSpace(mesh, elem, deg)
    sig = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
        
    # von Mises stress
    von_Mises = sqrt(3./2*inner(sig, sig))
    W = FunctionSpace(mesh, elem, deg)
    von_Mises = project(von_Mises, W)
        
    # Save to file and plot solution
    ufile << (u, t)
    mfile << (von_Mises, t)

 

Community: FEniCS Project
You use something like this:
cvtk = Function(V)
ph_f = File('./concentration.pvd')
ph_f << cvtk, t

# time loop
while t <= t_max:
     
       # assign variable
       cvtk.assign(p)
       ph_f << cvtk, t​
written 12 months ago by hirshikesh  
Thanks, but unfortunately that doesn't seem to work for me. If I say 
vM = Function(V)​

and then do my stress calculation in the time loop

# Calculate stress
sig = sigma(u) - (1./3)*tr(sigma(u))*Identity(d)  # deviatoric stress
      
# von Mises stress
von_Mises = sqrt(3./2*inner(sig, sig))
vM.assign(von_Mises)

I get the following error:

RuntimeError: 

*** Error:   Unable to function assignment.
*** Reason:  Expects a Function or linear combinations of Functions in the same FunctionSpaces.
*** Where:   This error was encountered inside function.py.
*** Process: 0
*** 
*** DOLFIN version: 2017.1.0
*** Git changeset:  70e77054c0a013e3a47e7161fa01ca1b83930c26
*** -------------------------------------------------------------------------
written 12 months ago by Alexandra K. Diem  
Why do you create a TensorFunctionSpace W which you never use? Also you don't need to create the FunctionSpace W in each time step.. Btw, you use W for three different things..
written 12 months ago by Adam Janecka  
Apologies, the TensorFunctionSpace was left over from when I tried something different and wanted to get the individual stress tensor components, I just forgot to remove it. Thanks for pointing out that I actually used W for both a Function and a variable, I absolutely did not notice that. I've been staring at this code for too long. It's the first time I'm trying to use Fenics beyond the tutorial problems, so please bear with me. 
written 12 months ago by Alexandra K. Diem  
No problem, everyone had to go through this ;-)
written 12 months ago by Adam Janecka  
can use:
z = FunctionSpace(mesh, 'CG',1)
cvtk = Function(z)​
written 11 months ago by hirshikesh  

1 Answer


6
12 months ago by
Before saving to the file, try to rename the von_Mises function, i.e.,
von_Mises.rename('vM', 'vM')​
Thanks, this solved the issue of having to rename the colour palettes after the simulation, but the value range in Paraview still varies for every time step.
written 12 months ago by Alexandra K. Diem  
In each time step, you can rescale the range to cover all the values ("Rescale to Data Range" button) -- this can be, of course, different for every time step. You can click on "Edit Color Map" and then you can either choose the data range manually (uncheck "Automatically Rescale to Fit Data Range" then "Rescale Range") or let ParaView rescale it over all time steps ("Rescale to Temporal Range").

Another option would be to normalize the values already in FEniCS which isn't probably the best option for a physical quantity like a von Mises stress.
written 12 months ago by Adam Janecka  
That works thank you! I was just wondering why this extra step is necessary for stress while it isn't for the solution variable on its own.
written 12 months ago by Alexandra K. Diem  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »