How to create correct VTK file for the material property


114
views
0
10 weeks 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)


F = u*v*dx + k*dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
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

1 Answer


2
10 weeks ago by
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/
Please login to add an answer/comment or follow this question.

Similar posts:
Search »