[PYTHON] Meshing polyhedra for micromagnetic simulation

background

For the simple micromagnetic simulation software MERRILL [^ 1] developed and released by Edinburgh University and others, we have summarized how to make a convex polyhedron into a simple mesh. The explanation of the mesh for MERRILL [^ 2] , </ sup> [^ 3] is limited to a very simple shape such as a rectangular parallelepiped, or an image data obtained by three-dimensionally scanning an actual sample. It seems that there is only English. With the method described here, you can create various shapes (such as the octahedron, which is the automorphic shape of magnetite) with CAD without using 3D scanning. I hope that students (and researchers) in geophysics and materials science will feel free to simulate and deepen their interest in magnetism. For those who are serious about it, it may be better to make a more decent mesh (I'm not sure).

0. Preparation

The compiled software of MERRILL is available on the Official Page. Use gmsh (4.7) to generate the mesh. I mainly used the api for python here. You can take it normally with pip. It's also in anaconda, but it doesn't seem to install properly ... The rest is file format conversion. MERRILL reads a mesh written in patran format or something, but I can't find an easy way to convert it. The community has created a matlab script called merrillsave, which is available ~~ published ~~ You can get it by emailing the author of the paper [^ 4]. It's a good idea to put it just for file format conversion, but I use octave to make it work. And it seems that this output is not exactly in patran format, and the shell script Convert2Pat.sh for converting again has been released [^ 2]. For file format conversion only (omitted), it is better to include wsl on Windows. The rest isn't important, but I used python (3.7) for some reason. Use meshio to extract the mesh information from the gmsh format file, and use scipy to convert it to a mat file. It is also convenient to remove the convex hull. paraview was used to visualize the calculation results.

1. Model generation

Prepare a polyhedral stl file. For example, let's make an octahedron using SketchUp. I made it as 200 (mm) on a side, but since the unit is arbitrary, I will treat it as 200 nm below.

octahedra-200nm.png

2. Mesh generation

It seems that gmsh has a mesh made from stl, so it may be easy to do. I will do it seriously here. In gmsh, points, lines connecting points, loops and faces created by lines, areas and volumes surrounded by faces are prepared in order, and finally, by applying mesh generation to the volume, faces and lines are saved. Make a mesh. You can use the points written in stl as they are, but since stl does not have line information itself, you have to make lines in reverse from the surface. Also, if you make a model properly, you often create a surface inside the polyhedron. Therefore, by removing the convex hull, only the outer surface is extracted. With these in mind, you could write, for example:

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

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

gmsh.initialize()

lc = 9 #Mesh size(9nm)

points = [gmsh.model.geo.addPoint(p[0], p[1], p[2], lc) for p in stldata.points] #Added points to gmsh. points is an array of point names

hull = ConvexHull(stldata.points) #Take the convex hull
simplices = np.array(hull.simplices)+1 #The index of each side of the convex hull. gmsh is 1-Since it is an index, add 1

spam = [] #Make a line from the surface. Omit the lines that are shared by the two sides.
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) #List of line endpoint pairs

lines = [gmsh.model.geo.addLine(l[0], l[1]) for l in edgepairs] #Add a line. Points are indicated by name. lines is an array of line names. You may wear dots and names.

def findcurve(l, edgepairs, lines): 
    '''
A function to define a face again with the extracted line.
Given the line l (end point pair), refer to the list of end point pairs and express it as the exponent of the saved line.
    '''
    b = (edgepairs[:, 0]==l[0]) & (edgepairs[:, 1]==l[1])
    if np.sum(b)>1:
        return 0 #There is no such thing as an index of zero, so this should result in an error somewhere.
    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) #The direction of the line is represented by a code.

loops = [] #There is a concept of loop in front of the face.
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] #Since it is a convex polyhedron, it defines faces for all loops.

gmsh.model.geo.addSurfaceLoop(surfaces, 100) #Since it is a polyhedron, all faces surround one volume. 100 is the name of the face loop (optional).
gmsh.model.geo.addVolume([100], 101) #Correspond the volume to the loop of the surface. 101 is the name of the volume (optional)

gmsh.model.geo.synchronize()
gmsh.model.mesh.generate(3) #3D mesh generation

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

Now you have a mesh. octahedra-200nm_3dmesh.png

(Display uses gmsh application)

Finally, leave the mesh in mat format for octave reading. Here, the points are stored as node, the tetrahedral mesh is stored as elem, and the triangular mesh is stored as face (although not used). Also, in MERRILL, the length unit is $ \ mathrm {\ mu} $ m.

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

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

#Since the matlab array starts from 1, add 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. File conversion

First on Octave

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

And then in Bash

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

Then, a pat file (like) was created.

4. Micromagnetic calculation

Let's calculate it as a trial. The parameters look like this

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

Save this as script.txt. As a magnetite of 20C, the initial state is homogeneous magnetization in the [111] direction, and the equilibrium state is searched for when a magnetic field of 0.05 mT is applied in the [111] direction. I saw the explanation that if you do something just on the axis of symmetry in numerical calculation, you often feel sick, so I'm shifting it a little. Perform the calculation

merrill script.txt

It took about 2 minutes (1000 steps).

The calculation result can be displayed in paraview.

result.png

It's vortex. I think that the z direction is the opposite in detail, but I'm not sure.

5. Failure (bonus)

I've made some mistakes before reaching this method, so I'll list them for reference.

According to the explanation [^ 3] of the Cambridge University group, 3D scan is set to stl by paraview, read as stlread on octave, and volume mesh is made because it is iso2mesh. Then, is it possible to output SketchUp stl like this ... First of all, stlread cannot read it. I'm not sure why this is, but in any case the polyhedral stl is so squishy that it's not a surface shape, so this probably won't work.

Since iso2mesh seems to use a C ++ library called CGAL, I also tried to make a volume mesh from stl with pygalmesh, which is a partial python port. I can do it (example with a 100nm octahedron),

octahedra-100nm_pygalmesh.png

Edges are not reproduced. I have a feeling that it will not affect the magnetic calculation (maybe it is a more natural form in a sense), but I am a little dissatisfied. It seems that you can specify the lines and faces to be saved in the original CGAL, but since pygalmesh is under development, I did not know what happened to that.

If it is an octahedron, you can create a 3D image numerically, but even if you use it, it will look the same. Real-life 3D scans deal with curved shapes, so edge preservation may not be a problem.

Recommended Posts

Meshing polyhedra for micromagnetic simulation