Einführung

In diesem Abschnitt wird beschrieben, wie die Temperaturdaten des zweidimensionalen Gitters im ndarray-Format, die durch die Wärmediffusionssimulation erhalten wurden, als VTK-Daten ausgegeben werden.

Methode

Ausgabe im StructuredGrid-Format mithilfe der Python VTK-Bibliothek.

Umgebung

--Windows 10 home (WSL wird für die FEniCS-Berechnung verwendet)

Anaconda(python 3.7.6)
FEniCS 2018

Simulationsergebnis

In der Temperaturdiffusionssimulation wurde die folgende Wärmediffusionsgleichung nach der Finite-Elemente-Methode unter Verwendung von FEniCS der Python-Schnittstelle berechnet. $ \frac{\partial}{\partial t} u = \alpha\nabla ^2 u $ $ u $: Temperaturskalar $ \ alpha $: Wärmediffusionskoeffizient. Das Berechnungsergebnis ist in der folgenden Abbildung dargestellt. Abbildung: Berechnungsergebnis

Informationen zu FEniCS finden Sie im folgenden Artikel. https://qiita.com/matsxxx/items/b696d2fe730be1683d8d Der Code zum Erstellen des Simulationsergebnisses, das dieses Mal verwendet wird, lautet wie folgt.

import fenics as fe
import numpy as np

tol = 1e-14 #
Lx = 80 #Länge in x-Richtung
Ly = 100 #Länge in y-Richtung
Nx = 80 #x-Achsen-Zellennummer
Ny = 100 #Anzahl der Zellen der y-Achse
T = 10.0 #Zeit zum Berechnen
num_steps = 10 #Zeitschritt
dt = T/num_steps #Zeitschritt
alpha = 100.#Wärmediffusionskoeffizient

#Netzerstellung
mesh = fe.RectangleMesh(fe.Point(0, 0), fe.Point(Lx, Ly), Nx, Ny)
V = fe.FunctionSpace(mesh, "P", 1)#Begrenzter Elementraum

#Wird verwendet, um die Reihenfolge der Fenics zu sortieren.
vertex_to_dof_map = fe.vertex_to_dof_map(V)

#Randbedingungspositionierung
def boundary_right(x, on_boundary):
    #x[0]:x Koordinaten, x[1]:y-Koordinate
    return on_boundary and fe.near(x[0], Lx, tol)


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

#Festlegung erstklassiger Randbedingungen
bc_right = fe.DirichletBC(V, fe.Constant(100), boundary_right)#100 ℃ auf der Lx-Seite
bc_top = fe.DirichletBC(V, fe.Constant(80), boundary_top)#80 ℃ auf Ly Seite

#Versuch/Testen Sie die Funktionseinstellungen
u = fe.TrialFunction(V) #Variablen, die Sie lösen möchten
v = fe.TestFunction(V) #Gewichtsfunktion

#Anfangsbedingungen
u_0 = fe.Expression('0',degree=2)#Anfangstemperatur 0
u_n = fe.interpolate(u_0, V)

#Variante der thermischen Diffusionsgleichung
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#Zeit

#Durchführung einer instationären Berechnung
for n in range(num_steps):
    #Zeitaktualisierung
    t += dt
    
    #Löse die Gleichung
    fe.solve(a == L, u, [bc_right, bc_top])
    
    #Lösung aktualisieren
    u_n.assign(u)

#Holen Sie sich das Endergebnis mit ndarray.
u_np_temp = u.vector().get_local()
#T aus der Liste der Variablen in Fenics[0,0], T[0,1]・ ・ ・ T.[Nx-1,Ny-2],T[Nx-1,Ny-1]Sequenzkonvertierung
u_np = np.zeros_like(u_np_temp)
for i,v in enumerate(vertex_to_dof_map):
    u_np[i] = u_np_temp[v]
    
#Speichern Sie mit npy
np.save("solution", u_np)

#Da es möglich ist, VTK-Daten an Fenics auszugeben, wird dies normalerweise verwendet.
#vtk = fe.File("output.pvd")
#vtk << u

Ausgabe mit VTK-Bibliothek

Das Simulationsergebnis im ndarray-Format wird im StructuredGrid-Format der VTK-Bibliothek ausgegeben.

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

#Simulationsergebnisse lesen
u = np.load("solution.npy")

#Erstellung von Punktdaten
points = vtk.vtkPoints()

for iy in range(Ny+1):
    for ix in range(Nx+1):
        x = ix 
        y = iy
        points.InsertNextPoint(x,y,0)#Da es zweidimensional ist, setzen Sie die z-Koordinate auf 0
    
#Sollwertdaten in strukturiertem Grid
structured_mesh = vtk.vtkStructuredGrid()
structured_mesh.SetExtent(0, Nx, 0, Ny, 0, 0)
structured_mesh.SetPoints(points)

#Geben Sie die Temperaturdaten in das vtk-Array ein
temp_value = vtk.vtkDoubleArray()
temp_value.SetName("Temperature")
for t in u:
    temp_value.InsertNextValue(t)
#Hinzufügen von Temperaturdaten zu Punktdaten des strukturierten Gitters
structured_mesh.GetPointData().AddArray(temp_value)

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

Bonus: Ausgabe mit PyEVTK

Wenn Sie nur als VTK-Daten ausgeben möchten, ist die Verwendung von PyEVTK einfacher.

import numpy as np 
import evtk

Nx = 80
Ny = 100
u = np.load("solution.npy")
u = u.reshape(Ny+1, -1)#x,Beachten Sie die Anordnung von y
u = u[:,:,np.newaxis]

#Erstellen Sie Koordinaten mit 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)#Da es zweidimensional ist, ist die Z-Koordinate 0

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