introduction

Cette section décrit comment sortir les données de température du réseau bidimensionnel au format ndarray obtenues par la simulation de diffusion thermique sous forme de données VTK.

Méthode

Sortie au format StructuredGrid à l'aide de la bibliothèque Python VTK.

environnement

--Windows 10 home (WSL est utilisé pour le calcul FEniCS)

Anaconda(python 3.7.6)
FEniCS 2018

résultat de la simulation

Dans la simulation de diffusion de température, l'équation de diffusion thermique suivante a été calculée par la méthode des éléments finis en utilisant FEniCS de l'interface Python. $ \frac{\partial}{\partial t} u = \alpha\nabla ^2 u $ $ u $: scalaire de température $ \ alpha $: coefficient de diffusion thermique. Le résultat du calcul est présenté dans la figure suivante. Figure: Résultat du calcul

Pour FEniCS, reportez-vous à l'article suivant. https://qiita.com/matsxxx/items/b696d2fe730be1683d8d Le code pour créer le résultat de la simulation utilisé cette fois est le suivant.

import fenics as fe
import numpy as np

tol = 1e-14 #
Lx = 80 #longueur de la direction x
Ly = 100 #longueur direction y
Nx = 80 #numéro de cellule sur l'axe des x
Ny = 100 #Nombre de cellules de l'axe y
T = 10.0 #Temps de calculer
num_steps = 10 #Pas de temps
dt = T/num_steps #Pas de temps
alpha = 100.#Coefficient de diffusion thermique

#Création de maillage
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(Lx, Ly), Nx, Ny)
V = fe.FunctionSpace(mesh, "P", 1)#Espace d'élément limité

#Utilisé pour trier la séquence des fenics.
vertex_to_dof_map = fe.vertex_to_dof_map(V)

#Positionnement des conditions aux limites
def boundary_right(x, on_boundary):
    #x[0]:coordonnées x, x[1]:coordonnée y
    return on_boundary and fe.near(x[0], Lx, tol)


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

#Définition de conditions aux limites de premier ordre
bc_right = fe.DirichletBC(V, fe.Constant(100), boundary_right)#100 ℃ côté Lx
bc_top = fe.DirichletBC(V, fe.Constant(80), boundary_top)#80 ℃ côté Ly

#Procès/Paramètres de la fonction de test
u = fe.TrialFunction(V) #Variables que vous souhaitez résoudre
v = fe.TestFunction(V) #Fonction poids

#Conditions initiales
u_0 = fe.Expression('0',degree=2)#Température initiale 0
u_n = fe.interpolate(u_0, V)

#Equation de diffusion thermique variable
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#temps

#Exécution d'un calcul non stationnaire
for n in range(num_steps):
    #Mise à jour de l'heure
    t += dt
    
    #Résous l'équation
    fe.solve(a == L, u, [bc_right, bc_top])
    
    #Mettre à jour la solution
    u_n.assign(u)

#Obtenez le résultat final avec ndarray.
u_np_temp = u.vector().get_local()
#T de la liste des variables en fenics[0,0], T[0,1]・ ・ ・ T[Nx-1,Ny-2],T[Nx-1,Ny-1]Conversion de séquence
u_np = np.zeros_like(u_np_temp)
for i,v in enumerate(vertex_to_dof_map):
    u_np[i] = u_np_temp[v]
    
#Économisez avec npy
np.save("solution", u_np)

#Puisqu'il est possible de sortir des données VTK vers des fenêtres, cela est généralement utilisé.
#vtk = fe.File("output.pvd")
#vtk << u

Sortie à l'aide de la bibliothèque VTK

Le résultat de la simulation au format ndarray est sorti au format StructuredGrid de la bibliothèque VTK.

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

#Lire les résultats de la simulation
u = np.load("solution.npy")

#Création de données ponctuelles
points = vtk.vtkPoints()

for iy in range(Ny+1):
    for ix in range(Nx+1):
        x = ix 
        y = iy
        points.InsertNextPoint(x,y,0)#Comme il est bidimensionnel, définissez la coordonnée z sur 0
    
#Données de point de consigne dans StructuredGrid
structured_mesh = vtk.vtkStructuredGrid()
structured_mesh.SetExtent(0, Nx, 0, Ny, 0, 0)
structured_mesh.SetPoints(points)

#Données de température d'entrée dans la baie vtk
temp_value = vtk.vtkDoubleArray()
temp_value.SetName("Temperature")
for t in u:
    temp_value.InsertNextValue(t)
#Ajouter des données de température aux données ponctuelles de la grille structurée
structured_mesh.GetPointData().AddArray(temp_value)

#Sortie de données
writer = vtk.vtkXMLDataObjectWriter()
writer.SetFileName("temperature.vts")
writer.SetInputData(structured_mesh)
writer.Update()

Bonus: sortie avec PyEVTK

Si vous souhaitez simplement générer des données VTK, il est plus facile d'utiliser PyEVTK.

import numpy as np 
import evtk

Nx = 80
Ny = 100
u = np.load("solution.npy")
u = u.reshape(Ny+1, -1)#x,Notez la disposition de y
u = u[:,:,np.newaxis]

#Créer des coordonnées avec 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)#Puisqu'elle est bidimensionnelle, la coordonnée Z est 0

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