### How to create correct VTK file for the material property

330
views
0
10 months ago by

I am trying to change the material property as time goes on and visualize it by Paraview. The code works; however, when I want to see the k.pvd in Paraview. It just shows the first step and for the next steps it gives me this error:

paraview: /build/paraview-arsa8T/paraview-5.0.1+dfsg1/VTK/IO/XML/vtkXMLDataReader.cxx:302: virtual void vtkXMLDataReader::SetupOutputData(): Assertion this->NumberOfPointArrays == this->PointDataArraySelection->GetNumberOfArraysEnabled()' failed.
Aborted (core dumped)


The code:

from __future__ import print_function
from fenics import *
import time

T = 2.0
num_steps = 5
dt = T / num_steps

# Create mesh and define function space
nx = ny = 10
mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)
V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition
def boundary(x, on_boundary):
return on_boundary

bc = DirichletBC(V, Constant(0), boundary)

# Define initial value
u_0 = Constant(10)
u_n = interpolate(u_0, V)

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Constant(0)
k1 = 1
k2 = 2
ref = 0
tol = 1e-7
k = Expression('x[0]> ref+tol ? k1 : k2',degree=0, k1=k1,k2=k2,tol=tol,ref=ref)

a, L = lhs(F), rhs(F)

vtkfile_T = File('heat_gaussian/T.pvd')
vtkfile_K = File('heat_gaussian/k.pvd')
u = Function(V)
t = 0.0
for n in range(num_steps):

t += dt
k.ref +=0.1
solve(a == L, u, bc)
vtkfile_T << (u, t)
vtkfile_K << (interpolate(k,V) ,t)
u_n.assign(u)

Community: FEniCS Project

I've got no error running your code. Probably, it depends on how your ParaView is compiled. Try to rename`, see https://bitbucket.org/fenics-project/dolfin/issues/80/reading-different-functions-with-different and https://www.allanswered.com/post/klnqb/saving-time-varying-expressions/