[PYTHON] Vernetzen Sie das Polyeder für die mikromagnetische Simulation

Hintergrund

Für die einfache mikromagnetische Simulationssoftware MERRILL [^ 1], die von der Universität Edinburgh und anderen entwickelt und veröffentlicht wurde, haben wir zusammengefasst, wie aus einem konvexen Polyeder ein einfaches Netz gemacht werden kann. Die Erklärung des Netzes für MERRILL [^ 2] , </ sup> [^ 3] beschränkt sich auf eine sehr einfache Form wie einen rechteckigen Körper oder eine Form, die Bilddaten verwendet, die durch dreidimensionales Scannen einer tatsächlichen Probe erhalten wurden. Es scheint, dass es nur Englisch gibt. Ich denke, dass Sie mit der hier beschriebenen Methode verschiedene Formen (wie Oktaeder, ein selbstgeformtes magnetisches Eisenerz) mit CAD ohne dreidimensionales Scannen herstellen können. Ich hoffe, dass Studenten (und Forscher) in Geophysik und Materialwissenschaften ihr Interesse am Magnetismus simulieren und vertiefen können. Für diejenigen, die es ernst meinen, ist es möglicherweise besser, ein anständigeres Netz herzustellen (ich bin mir nicht sicher).

0. Vorbereitung

Die kompilierte Software von MERRILL ist auf der [offiziellen Seite] verfügbar (https://www.geos.ed.ac.uk/geosciences/research/projects/rockmag). Verwenden Sie gmsh (4.7), um das Netz zu generieren. Ich habe hier hauptsächlich API für Python verwendet. Sie können es normal mit pip nehmen. Es ist auch in Anaconda, aber es scheint nicht richtig zu installieren ... Der Rest ist die Konvertierung des Dateiformats. MERRILL liest Netze, die im Patran-Format usw. geschrieben sind, aber ich kann keine einfache Möglichkeit finden, sie zu konvertieren. Ein Matlab-Skript namens merrill save wurde von der Community erstellt und kann per E-Mail an den Autor des Artikels [^ 4] bezogen werden. Es ist eine gute Idee, es nur für die Dateiformatkonvertierung zu verwenden, aber ich verwende Oktave, damit es funktioniert. Und es scheint, dass diese Ausgabe nicht genau im Patran-Format vorliegt und das Shell-Skript Convert2Pat.sh zum erneuten Konvertieren veröffentlicht wurde [^ 2]. Nur für Windows sollte wsl nur für die Dateiformatkonvertierung enthalten sein (weggelassen). Der Rest ist nicht wichtig, aber ich habe aus irgendeinem Grund Python (3.7) verwendet. Verwenden Sie meshio, um die Netzinformationen aus der Datei im gmsh-Format zu extrahieren, und konvertieren Sie sie mit scipy in eine Mat-Datei. Es ist auch bequem, ein konvexes Paket zu nehmen. paraview wurde verwendet, um die Berechnungsergebnisse zu visualisieren.

1. Modellgenerierung

Bereiten Sie eine Polyeder-STL-Datei vor. Lassen Sie uns zum Beispiel mit SketchUp ein Oktaeder erstellen. Ich habe es auf einer Seite als 200 (mm) gemacht, aber da die Einheit willkürlich ist, werde ich es unten als 200 nm behandeln.

octahedra-200nm.png

2. Netzgenerierung

Es scheint, dass gmsh ein Netz aus stl hat, so dass es einfach sein kann. Ich werde es hier ernsthaft tun. In gmsh werden Punkte, Linien, die Punkte, Schleifen und Flächen verbinden, die durch Linien, Bereiche und Volumina erzeugt werden, die von Flächen umgeben sind, in der richtigen Reihenfolge vorbereitet. Schließlich wird die Netzgenerierung auf das Volume angewendet, um Flächen und Linien zu speichern. Machen Sie ein Netz. Sie können die in stl geschriebenen Punkte so verwenden, wie sie sind, aber stl verfügt nicht über die Linieninformationen selbst, sodass Sie eine Linie in umgekehrter Reihenfolge von der Oberfläche erstellen müssen. Wenn Sie ein Modell richtig erstellen, erstellen Sie häufig eine Oberfläche innerhalb des Polyeders. Daher wird durch Entfernen des konvexen Pakets nur die äußere Oberfläche extrahiert. In diesem Sinne könnten Sie zum Beispiel schreiben:

import gmsh
import meshio
import numpy as np
from scipy.spatial import ConvexHull

stldata = meshio.read("octahedra-200nm.stl")

gmsh.initialize()

lc = 9 #Maschenweite(9nm)

points = [gmsh.model.geo.addPoint(p[0], p[1], p[2], lc) for p in stldata.points] #Punkte zu gmsh hinzugefügt. Punkte ist ein Array von Punktnamen

hull = ConvexHull(stldata.points) #Nehmen Sie ein konvexes Paket
simplices = np.array(hull.simplices)+1 #Der Index jeder Seite des konvexen Pakets. gmsh ist 1-Da es sich um einen Index handelt, addieren Sie 1

spam = [] #Machen Sie eine Linie von der Oberfläche. Lassen Sie die Zeilen weg, die von beiden Seiten geteilt werden.
for l in simplices:
    spam.extend([list({l[0], l[1]}), list({l[1], l[2]}), list({l[2], l[0]})])
edgepairs = []
for l in spam:
    if l not in edgepairs:
        edgepairs.append(l)
edgepairs = np.array(edgepairs) #Liste der Linienendpunktpaare

lines = [gmsh.model.geo.addLine(l[0], l[1]) for l in edgepairs] #Fügen Sie eine Zeile hinzu. Punkte werden namentlich angegeben. Zeilen ist ein Array von Zeilennamen. Punkte und Namen können getragen werden.

def findcurve(l, edgepairs, lines): 
    '''
Eine Funktion zum erneuten Definieren eines Gesichts mit den extrahierten Linien.
Beziehen Sie sich angesichts der Linie l (Endpunktpaar) auf die Liste der Endpunktpaare und drücken Sie sie als Exponenten der gespeicherten Linie aus.
    '''
    b = (edgepairs[:, 0]==l[0]) & (edgepairs[:, 1]==l[1])
    if np.sum(b)>1:
        return 0 #Es gibt keinen Index von Null, daher sollte dies irgendwo zu einem Fehler führen.
    elif np.sum(b)==1:
        return np.where(b==True)[0][0]+1
    else:
        b = (edgepairs[:, 0]==l[1]) & (edgepairs[:, 1]==l[0])
        if np.sum(b)==0:
            return 0
        else:
            #print(np.where(b==True)[0][0])
            return -1 * (np.where(b==True)[0][0]+1) #Die Richtung der Linie wird durch einen Code dargestellt.

loops = [] #Es gibt ein Konzept der Schleife vor dem Gesicht.
for s in simplices:
    l1 = findcurve([s[0], s[1]], edgepairs, lines)
    l2 = findcurve([s[1], s[2]], edgepairs, lines)
    l3 = findcurve([s[2], s[0]], edgepairs, lines)
    print([l1, l2, l3])
    loops.append(gmsh.model.geo.addCurveLoop([l1, l2, l3]))

surfaces = [gmsh.model.geo.addPlaneSurface([loop]) for loop in loops] #Da es sich um ein konvexes Polyeder handelt, werden Flächen für alle Schleifen definiert.

gmsh.model.geo.addSurfaceLoop(surfaces, 100) #Da es sich um ein Polyeder handelt, umgeben alle Flächen ein Volumen. 100 ist der Name der Gesichtsschleife (optional).
gmsh.model.geo.addVolume([100], 101) #Entsprechen Sie dem Volumen der Schleife der Oberfläche. 101 ist der Name des Volumes (optional)

gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3) #Erzeugung eines 3D-Netzes

gmsh.write("octahedra-200nm.msh")
gmsh.finalize()

Jetzt hast du ein Netz. octahedra-200nm_3dmesh.png

(Anzeige verwendet gmsh-Anwendung)

Lassen Sie das Netz schließlich im Mattenformat, um die Oktave zu lesen. Hier werden die Punkte als Knoten gespeichert, das Tetraedernetz wird als Element gespeichert und das Dreiecksnetz wird als Fläche gespeichert (obwohl nicht verwendet). In MERRILL ist die Längeneinheit $ \ mathrm {\ mu} $ m.

import meshio
from scipy.io import savemat
m = meshio.read("octahedra-200nm.msh")

# nm -> um
node = m.points / 1000.

#Da das Matlab-Array bei 1 beginnt, addieren Sie 1
elem = m.cells_dict["tetra"] + 1
face = m.cells_dict["triangle"] + 1
spam = {"node": node, "elem": elem, "face": face}

savemat("octahedra-200nm.mat", spam)

3. Dateikonvertierung

Zuerst auf Oktave

load("octahedra-200nm.mat")
merrillsavepat("octahedra-200nm_K", node, elem)

Und dann in Bash

Convert2Pat.sh octahedra-200nm_K octahedra-200nm.pat

Dann wurde eine Pat-Datei (wie) erstellt.

4. Mikromagnetische Berechnung

Berechnen wir es als Versuch. Die Parameter sehen so aus

Set MaxMeshNumber 1

ReadMesh 1 octahedra-200nm.pat

Set MaxEnergyEvaluations 5000

ConjugateGradient
Set ExchangeCalculator 1

Magnetite 20 C
Uniform magnetization 0.99,1,1
External Field Strength 0.05 mT
External Field Direction 1,1,0.999

EnergyLog energy_log.txt

Minimize
WriteMagnetization Octahedra-200nm_log
CloseLogfile
END

Speichern Sie dies als script.txt. Als magnetisches Eisenerz mit 20 ° C ist der Anfangszustand eine gleichmäßige Magnetisierung in [111] -Richtung, und der Gleichgewichtszustand wird gesucht, wenn ein Magnetfeld von 0,05 mT in [111] -Richtung angelegt wird. Ich habe die Erklärung gesehen, dass man sich oft krank fühlt, wenn man etwas nur auf der Symmetrieachse numerisch macht, also verschiebe ich es ein wenig. Führen Sie die Berechnung durch

merrill script.txt

Es dauerte ungefähr 2 Minuten (1000 Schritte).

Das Berechnungsergebnis kann in der Paraview angezeigt werden.

result.png

Es ist Wirbel. Ich denke, dass die z-Richtung im Detail das Gegenteil ist, aber ich bin mir nicht sicher.

5. Misserfolg (Bonus)

Ich habe einige Fehler gemacht, bevor ich diese Methode erreicht habe, daher werde ich sie als Referenz auflisten.

Gemäß der Erklärung [^ 3] der großen Gruppe von Cambridge wird der dreidimensionale Scan durch Paraview in stl umgewandelt, als stlread auf Oktave gelesen und als Volumennetz erstellt, da er als iso2mesh bezeichnet wird. Kann die stl-Ausgabe von SketchUp dann so ausgeführt werden? Erstens kann sie möglicherweise nicht von stlread gelesen werden. Ich bin mir nicht sicher, warum das so ist, aber auf jeden Fall ist der stl des Polyeders so matschig, dass es nicht als Oberflächenform bezeichnet werden kann, also wird dies wahrscheinlich nicht funktionieren.

Da iso2mesh anscheinend eine C ++ - Bibliothek namens CGAL verwendet, habe ich auch versucht, aus stl mit pygalmesh, einem partiellen Python-Port, ein Volume-Mesh zu erstellen. Ich kann es tun (Beispiel mit 100nm Oktaeder),

octahedra-100nm_pygalmesh.png

Kanten werden nicht reproduziert. Ich habe das Gefühl, dass es die magnetische Berechnung nicht beeinflusst (vielleicht ist es in gewissem Sinne eine natürlichere Form), aber ich bin ein wenig unzufrieden. Es scheint, dass Sie die Linien und Flächen angeben können, die in der ursprünglichen CGAL gespeichert werden sollen, aber da sich Pygalmesh in der Entwicklung befindet, wusste ich nicht, was passiert ist.

Wenn es sich um ein Oktaeder handelt, können Sie ein dreidimensionales Bild numerisch erstellen, aber selbst wenn Sie es verwenden, sieht es gleich aus. Die Kantenerhaltung ist möglicherweise kein großes Problem, da reale dreidimensionale Scans sich mit gekrümmten Formen befassen.

Recommended Posts

Vernetzen Sie das Polyeder für die mikromagnetische Simulation