Introduction

This section describes how to output the temperature data of the 2D grid in ndarray format obtained by the thermal diffusion simulation as VTK data.

Method

Output in StructuredGrid format using Python VTK library.

environment

--Windows 10 home (WSL is used for FEniCS calculation)

Anaconda(python 3.7.6)
FEniCS 2018

simulation result

In the thermal diffusion simulation, the following thermal diffusion equation was calculated by the finite element method using FEniCS of the Python interface. $ \frac{\partial}{\partial t} u = \alpha\nabla ^2 u $ $ u $: Temperature scalar $ \ alpha $: Thermal diffusivity. The calculation result is shown in the following figure. Figure: Calculation result

For FEniCS, refer to the following article. https://qiita.com/matsxxx/items/b696d2fe730be1683d8d The code to create the simulation result used this time is as follows.

import fenics as fe
import numpy as np

tol = 1e-14 #
Lx = 80 #x-direction length
Ly = 100 #y-direction length
Nx = 80 #x-axis cell count
Ny = 100 #Number of y-axis cells
T = 10.0 #Time to calculate
num_steps = 10 #Time step
dt = T/num_steps #Time step
alpha = 100.#Thermal diffusion coefficient

#Mesh creation
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(Lx, Ly), Nx, Ny)
V = fe.FunctionSpace(mesh, "P", 1)#Finite element space

#Used to sort the fenics array.
vertex_to_dof_map = fe.vertex_to_dof_map(V)

#Boundary condition positioning
def boundary_right(x, on_boundary):
    #x[0]:x coordinate, x[1]:y coordinate
    return on_boundary and fe.near(x[0], Lx, tol)


def boundary_top(x, on_boundary):
    #x[0]:x coordinate, x[1]:y coordinate
    return on_boundary and fe.near(x[1], Ly, tol)

#Setting of first-class boundary conditions
bc_right = fe.DirichletBC(V, fe.Constant(100), boundary_right)#100 ℃ on the Lx side
bc_top = fe.DirichletBC(V, fe.Constant(80), boundary_top)#80 ℃ on Ly side

#Trial/Test function settings
u = fe.TrialFunction(V) #Variables you want to solve
v = fe.TestFunction(V) #Weight function

#Initial conditions
u_0 = fe.Expression('0',degree=2)#Initial temperature 0
u_n = fe.interpolate(u_0, V)

#Variational thermal diffusion equation
F= u * v * fe.dx + alpha * dt * fe.dot(fe.grad(u), fe.grad(v)) * fe.dx - (u_n) * v * fe.dx 
a, L = fe.lhs(F), fe.rhs(F)

#variable
u = fe.Function(V)
t = 0#time

#Execution of non-stationary calculation
for n in range(num_steps):
    #Time update
    t += dt
    
    #Solve the equation
    fe.solve(a == L, u, [bc_right, bc_top])
    
    #Update solution
    u_n.assign(u)

#Get the final result with ndarray.
u_np_temp = u.vector().get_local()
#T from the list of variables in fenics[0,0], T[0,1]・ ・ ・ T[Nx-1,Ny-2],T[Nx-1,Ny-1]Sequence conversion
u_np = np.zeros_like(u_np_temp)
for i,v in enumerate(vertex_to_dof_map):
    u_np[i] = u_np_temp[v]
    
#Save with npy
np.save("solution", u_np)

#Since it is possible to output VTK data to fenics, this is usually used.
#vtk = fe.File("output.pvd")
#vtk << u

Output using VTK library

The simulation result in ndarray format is output in StructuredGrid format of VTK library.

import numpy as np
import vtk
Nx = 80
Ny = 100

#Read simulation results
u = np.load("solution.npy")

#Creation of point data
points = vtk.vtkPoints()

for iy in range(Ny+1):
    for ix in range(Nx+1):
        x = ix 
        y = iy
        points.InsertNextPoint(x,y,0)#Since it is two-dimensional, the z coordinate should be 0.
    
#Set point data in structuredGrid
structured_mesh = vtk.vtkStructuredGrid()
structured_mesh.SetExtent(0, Nx, 0, Ny, 0, 0)
structured_mesh.SetPoints(points)

#Input temperature data to vtk array
temp_value = vtk.vtkDoubleArray()
temp_value.SetName("Temperature")
for t in u:
    temp_value.InsertNextValue(t)
#Add temperature data to Point Data of Structured Grid
structured_mesh.GetPointData().AddArray(temp_value)

#Data output
writer = vtk.vtkXMLDataObjectWriter()
writer.SetFileName("temperature.vts")
writer.SetInputData(structured_mesh)
writer.Update()

Bonus: Output using PyEVTK

If you just want to output as VTK data, it is easier to use PyEVTK.

import numpy as np 
import evtk

Nx = 80
Ny = 100
u = np.load("solution.npy")
u = u.reshape(Ny+1, -1)#x,Note the arrangement of y
u = u[:,:,np.newaxis]

#Create coordinates with meshgrid
x = np.arange(0, Nx+1)
y = np.arange(0, Ny+1)
X, Y = np.meshgrid(x,y)
X = X[:,:,np.newaxis]
Y = Y[:,:,np.newaxis]

Z = np.zeros_like(X)#Since it is two-dimensional, the Z coordinate is 0.

evtk.hl.gridToVTK("temperature_evtk", X, Y, Z, pointData={"temperature":u})