[PYTHON] Mailler le polyèdre pour la simulation micromagnétique

Contexte

Pour le logiciel de simulation micro magnétique simple MERRILL [^ 1] développé et publié par l'Université d'Édimbourg et d'autres, nous avons résumé comment transformer un polyèdre convexe en un maillage simple. L'explication du maillage pour MERRILL [^ 2] , </ sup> [^ 3] est limitée à une forme très simple comme un corps rectangulaire ou qui utilise des données d'image obtenues par balayage tridimensionnel d'un échantillon réel. Il semble qu'il n'y ait que l'anglais. Avec la méthode décrite ici, je pense que vous pouvez créer diverses formes (comme l'octaèdre, qui est un minerai de fer magnétique en forme de soi) avec CAD sans balayage tridimensionnel. J'espère que les étudiants (et chercheurs) en géophysique et science des matériaux se sentiront libres de simuler et d'approfondir leur intérêt pour le magnétisme. Pour ceux qui sont sérieux, il vaut peut-être mieux faire un maillage plus décent (je ne suis pas sûr).

0. Préparation

Le logiciel compilé de MERRILL est disponible sur la page officielle. Utilisez gmsh (4.7) pour générer le maillage. J'ai principalement utilisé l'api pour python ici. Vous pouvez le prendre normalement avec pip. C'est aussi dans anaconda, mais il ne semble pas s'installer correctement ... Le reste est la conversion de format de fichier. MERRILL lit les maillages écrits au format patran et ainsi de suite, mais je ne trouve pas de moyen facile de les convertir. Un script matlab appelé merrill save a été créé par la communauté et peut être obtenu en envoyant un email à l'auteur de l'article [^ 4]. C'est une bonne idée de le mettre juste pour la conversion de format de fichier, mais j'utilise octave pour le faire fonctionner. Et il semble que cette sortie ne soit pas exactement au format patran, et le script shell Convert2Pat.sh pour une nouvelle conversion a été publié [^ 2]. Pour Windows uniquement, wsl doit être inclus uniquement pour la conversion de format de fichier (omis). Le reste n'est pas important, mais j'ai utilisé python (3.7) pour une raison quelconque. Utilisez meshio pour extraire les informations de maillage du fichier au format gmsh et utilisez scipy pour les convertir en fichier mat. Il est également pratique de prendre un emballage convexe. paraview a été utilisé pour visualiser les résultats des calculs.

1. Génération de modèles

Préparez un fichier stl polyèdre. Par exemple, faisons un octaèdre en utilisant SketchUp. Je l'ai fait à 200 (mm) d'un côté, mais comme l'unité est arbitraire, je la traiterai comme 200 nm en dessous.

octahedra-200nm.png

2. Génération de maillage

Il semble que gmsh a un maillage fait de stl, donc cela peut être facile à faire. Je vais le faire sérieusement ici. Dans gmsh, les points, les lignes reliant les points, les boucles et les faces créés par des lignes, des zones et des volumes entourés de faces sont préparés dans l'ordre, et enfin, la génération de maillage est appliquée au volume pour enregistrer les faces et les lignes. Faites un maillage. Vous pouvez utiliser les points écrits en stl tels quels, mais stl n'a pas les informations de ligne elle-même, vous devez donc créer une ligne à l'envers de la surface. De plus, si vous créez correctement un modèle, il arrive souvent qu'une surface soit créée à l'intérieur du polyèdre. Par conséquent, en retirant l'emballage convexe, seule la surface extérieure est extraite. Dans cet esprit, vous pouvez écrire, par exemple:

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

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

gmsh.initialize()

lc = 9 #Mesh taille(9nm)

points = [gmsh.model.geo.addPoint(p[0], p[1], p[2], lc) for p in stldata.points] #Ajout de points à gmsh. points est un tableau de noms de points

hull = ConvexHull(stldata.points) #Prenez un emballage convexe
simplices = np.array(hull.simplices)+1 #L'index de chaque côté du paquet convexe. gmsh vaut 1-Puisqu'il s'agit d'un index, ajoutez 1

spam = [] #Tracez une ligne à partir de la surface. Omettez les lignes partagées par les deux côtés.
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 des paires de points d'extrémité de ligne

lines = [gmsh.model.geo.addLine(l[0], l[1]) for l in edgepairs] #Ajoutez une ligne. Les points sont indiqués par leur nom. lines est un tableau de noms de lignes. Des points et des noms peuvent être portés.

def findcurve(l, edgepairs, lines): 
    '''
Une fonction pour redéfinir un visage avec les lignes extraites.
Étant donné la ligne l (paire de points d'extrémité), reportez-vous à la liste des paires de points d'extrémité et exprimez-la comme l'exposant de la ligne enregistrée.
    '''
    b = (edgepairs[:, 0]==l[0]) & (edgepairs[:, 1]==l[1])
    if np.sum(b)>1:
        return 0 #Il n'existe pas d'index de zéro, cela devrait donc entraîner une erreur quelque part.
    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) #La direction de la ligne est représentée par un code.

loops = [] #Il existe un concept de boucle devant le visage.
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] #Puisqu'il s'agit d'un polyèdre convexe, les faces sont définies pour toutes les boucles.

gmsh.model.geo.addSurfaceLoop(surfaces, 100) #Puisqu'il s'agit d'un polyèdre, toutes les faces entourent un volume. 100 est le nom de la boucle de visage (facultatif).
gmsh.model.geo.addVolume([100], 101) #Faites correspondre le volume à la boucle de la surface. 101 est le nom du volume (facultatif)

gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3) #Génération de maillage 3D

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

Vous avez maintenant un maillage. octahedra-200nm_3dmesh.png

(L'affichage utilise l'application gmsh)

Enfin, laissez le maillage au format mat pour la lecture d'octave. Ici, les points sont stockés en tant que nœud, le maillage tétraédrique est stocké en tant qu'élem et le maillage triangulaire est stocké en tant que face (bien que non utilisé). De plus, dans MERRILL, l'unité de longueur est $ \ mathrm {\ mu} $ m.

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

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

#Puisque le tableau matlab commence à 1, ajoutez 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. Conversion de fichiers

Premier sur Octave

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

Et puis dans Bash

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

Ensuite, un fichier pat (comme) a été créé.

4. Calcul micro magnétique

Calculons-le comme un essai. Les paramètres ressemblent à ceci

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

Enregistrez-le sous script.txt. En tant que minerai de fer magnétique 20C, l'état initial est une aimantation uniforme dans la direction [111], et l'état d'équilibre est recherché lorsqu'un champ magnétique de 0,05 mT est appliqué dans la direction [111]. J'ai vu l'explication que si vous faites quelque chose juste sur l'axe de symétrie numériquement, vous vous sentez souvent malade, alors je change un peu. Effectuer le calcul

merrill script.txt

Cela a pris environ 2 minutes (1000 étapes).

Le résultat du calcul peut être affiché dans paraview.

result.png

C'est un vortex. Je pense que la direction z est le contraire en détail, mais je ne suis pas sûr.

5. Échec (bonus)

J'ai fait quelques erreurs avant d'arriver à cette méthode, je vais donc les énumérer pour référence.

Selon le commentaire [^ 3] du Cambridge Large Group, le scan tridimensionnel est stl en paraview, lu comme stlread sur octave, et fait un maillage de volume car il est iso2mesh. Ensuite, est-il possible de sortir la sortie stl de SketchUp comme ça? Je ne sais pas pourquoi, mais dans tous les cas, la stl polyédrique est si molle qu'elle n'est pas une forme de surface, donc cela ne fonctionnera probablement pas.

Comme iso2mesh semble utiliser une bibliothèque C ++ appelée CGAL, j'ai également essayé de créer un maillage de volume à partir de stl avec pygalmesh, qui est un portage python partiel. Je peux le faire (exemple avec l'octaèdre 100nm),

octahedra-100nm_pygalmesh.png

Les bords ne sont pas reproduits. J'ai le sentiment que cela n'affectera pas le calcul magnétique (c'est peut-être une forme plus naturelle dans un sens), mais je suis un peu insatisfait. Il semble que vous puissiez spécifier les lignes et les faces à enregistrer dans le CGAL d'origine, mais comme pygalmesh est en cours de développement, je ne savais pas ce qui s'était passé.

Vous pouvez créer une image tridimensionnelle numériquement s'il s'agit d'un octaèdre, mais même si vous l'utilisez, elle aura le même aspect. La préservation des bords n'est peut-être pas un problème majeur, car les scans tridimensionnels réels traitent des formes courbes.

Recommended Posts

Mailler le polyèdre pour la simulation micromagnétique