Feature: GeometryEdgeFinder

doc/templates/api-template.html

@@ -42,7 +42,7 @@ <h3>Basics</h3>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'mesh')">Mesh size and geometry</a></li>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'shapes')">Shapes</a></li>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'regions')">Material regions</a></li>
-    <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'initial')">Initial Magnetization</a></li>
+    <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'initial')">Initial magnetization</a></li>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'parameters')">Material parameters</a></li>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'excitation')">Excitation</a></li>
     <li class="pure-menu-item"><a href="#" class="pure-menu-link api-section-link" onclick="openSection(event, 'outputquantities')">Output quantities</a></li>

doc/templates/examples-template.html

@@ -222,7 +222,7 @@ <h2 id="ex_SP4">Getting started with Standard Problem #4</h2>
 
 Note: these are 3D geometries seen from above. The displayed cell filling is averaged along the thickness (notable in ellipse and layers example). Black means empty space, white is filled.
 
-<hr/><h2 id="ex_initmag">Initial Magnetization</h2>
+<hr/><h2 id="ex_initmag">Initial magnetization</h2>
 Some initial magnetization functions are provided, as well as transformations similar to those on Shapes. See the Config <a href="http://mumax.github.io/api.html">API</a>.
 
 {{.Example `
@@ -280,7 +280,7 @@ <h2 id="ex_SP4">Getting started with Standard Problem #4</h2>
 The magnetization can also be set in separate regions, see below.
 
 
-<hr/><h2 id="ex_cheese">Interlude: Rotating Cheese</h2>
+<hr/><h2 id="ex_cheese">Interlude: rotating cheese</h2>
 
 In this example we define a geometry that looks like a slice of cheese and have it rotate in time.
 
@@ -317,7 +317,7 @@ <h2 id="ex_SP4">Getting started with Standard Problem #4</h2>
 {{.Output}}
 
 
-<hr/><h2 id="ex_regions">Regions: Space-dependent Parameters</h2>
+<hr/><h2 id="ex_regions">Regions: space-dependent parameters</h2>
 
 <p>Space-dependent parameters are defined using material <i>regions</i>. Regions are numbered 0-255 and represent different materials. Each cell can belong to only one region. At the start of a simulation all cells have region number 0.</p>
 
@@ -429,7 +429,7 @@ <h2 id="ex_SP4">Getting started with Standard Problem #4</h2>
 {{.Output}}
 
 
-<hr/><h2 id="ex_PMA">PMA Racetrack</h2>
+<hr/><h2 id="ex_PMA">PMA racetrack</h2>
 In this example we drive a domain wall in PMA material by spin-transfer torque.  We set up a post-step function that makes the simulation box "follow" the domain wall. Like this, only a small number of cells is needed to simulate an infinitely long magnetic wire.
 
 {{.Example `
@@ -463,7 +463,7 @@ <h2 id="ex_SP4">Getting started with Standard Problem #4</h2>
 Since we center on the domain wall we can not see that it is actually moving, but the domain wall breakdown is visible.
 
 
-<hr/><h2 id="ex_Py">Py Racetrack</h2>
+<hr/><h2 id="ex_Py">Py racetrack</h2>
 
 In this example we drive a vortex wall in Permalloy by spin-transfer torque. The simulation box "follows" the domain wall. By removing surface charges at the left and right ends, we mimic an infintely long wire.
 

engine/ext_geomedgefinder.go

@@ -0,0 +1,102 @@
+package engine
+
+import (
+	"math"
+	"strings"
+
+	"github.com/mumax/3/cuda"
+	"github.com/mumax/3/util"
+)
+
+func init() {
+	DeclFunc("ext_GeomEdge", ext_GeomEdge, "Returns an edge (+x, -x, +y, -y, +z, -z) of the geometry as a Shape. The Shape may not update automatically after changes to the mesh or geometry; rerun the function to refresh.")
+}
+
+func ext_GeomEdge(axis string) Shape {
+
+	s, r := geometry.Slice()
+	if r {
+		defer cuda.Recycle(s)
+	}
+
+	arr3d := s.HostCopy().Scalars()
+
+	n := Mesh().Size()
+	Nx, Ny, Nz := n[X], n[Y], n[Z]
+	PBC := Mesh().PBC()
+
+	edgemask := make([]bool, Nx*Ny*Nz)
+
+	var offx, offy, offz int
+	switch strings.ToLower(axis) {
+	case "+x":
+		offx, offy, offz = 1, 0, 0
+	case "-x":
+		offx, offy, offz = -1, 0, 0
+	case "+y":
+		offx, offy, offz = 0, 1, 0
+	case "-y":
+		offx, offy, offz = 0, -1, 0
+	case "+z":
+		offx, offy, offz = 0, 0, 1
+	case "-z":
+		offx, offy, offz = 0, 0, -1
+	default:
+		util.Log("GeometryEdge: invalid direction '" + axis + "'. Use +x, -x, +y, -y, +z, or -z")
+	}
+
+	for k := 0; k < Nz; k++ {
+		for j := 0; j < Ny; j++ {
+			for i := 0; i < Nx; i++ {
+
+				if arr3d[k][j][i] == 0 {
+					continue
+				}
+
+				wrapPBC := func(i, N, PBC int) int {
+					if PBC == 0 {
+						return i // Don't wrap
+					}
+					return ((i % N) + N) % N // Wrap in integer range [0, N-1]
+				}
+
+				nx, ny, nz := wrapPBC(i+offx, Nx, PBC[X]), wrapPBC(j+offy, Ny, PBC[Y]), wrapPBC(k+offz, Nz, PBC[Z])
+				if nx < 0 || nx >= Nx || ny < 0 || ny >= Ny || nz < 0 || nz >= Nz || arr3d[nz][ny][nx] == 0 {
+					edgemask[(k*Ny+j)*Nx+i] = true
+				}
+			}
+		}
+	}
+
+	if PBC[X] != 0 || PBC[Y] != 0 || PBC[Z] != 0 { // Account for PBC in case of Grid resize by repeating GeomEdge
+		d := Mesh().CellSize()
+		Rx := float64(Nx) * d[X] * float64(PBC[X]&1)
+		Ry := float64(Ny) * d[Y] * float64(PBC[Y]&1)
+		Rz := float64(Nz) * d[Z] * float64(PBC[Z]&1)
+		return maskToShape(edgemask).Repeat(Rx, Ry, Rz)
+	} else { // No PBC
+		return maskToShape(edgemask)
+	}
+}
+
+// Helper function to wrap bool mask into a Shape
+func maskToShape(mask []bool) Shape {
+	n := Mesh().Size()
+	Nx, Ny, Nz := n[X], n[Y], n[Z]
+	d := Mesh().CellSize()
+
+	return func(x, y, z float64) bool {
+		Lx := float64(Nx) * d[X]
+		Ly := float64(Ny) * d[Y]
+		Lz := float64(Nz) * d[Z]
+
+		ix := int(math.Floor((x + 0.5*Lx) / d[X]))
+		iy := int(math.Floor((y + 0.5*Ly) / d[Y]))
+		iz := int(math.Floor((z + 0.5*Lz) / d[Z]))
+
+		if ix < 0 || ix >= Nx || iy < 0 || iy >= Ny || iz < 0 || iz >= Nz {
+			return false
+		}
+		return mask[(iz*Ny+iy)*Nx+ix]
+	}
+}

test/ext_GeometryEdge.mx3

@@ -0,0 +1,72 @@
+/*
+    This test primarily serves as a visual verification of the ext_geomEdge() functionality,
+    but also performs two automated tests for particular analytically describable cases with PBC along X and Z.
+*/
+
+N := 21
+c := 5e-9
+d := 115e-9
+
+SetCellSize(c, c, c)
+SetGridSize(N, N, N)
+SetPBC(1,0,1) // Mix of PBC and OBC to simultaneously check behavior for both cases
+
+m = Uniform(-1, 0, 0)
+shape := SuperBall(d, 1.5) // Shape with analytical expression available
+SetGeom(shape)
+SaveAs(m, "m_orig.ovf")
+
+PlusX := ext_GeomEdge("+X")
+PlusY := ext_GeomEdge("+Y")
+MinZ := ext_GeomEdge("-Z")
+
+SetGeom(PlusX)
+SaveAs(m, "m_+X.ovf")
+SetGeom(PlusY)
+SaveAs(m, "m_+Y.ovf")
+SetGeom(MinZ)
+SaveAs(m, "m_-Z.ovf")
+
+//// Double grid size to see effect of PBC
+SetGridSize(2*N+1, 2*N+1, 2*N+1)
+
+SetGeom(PlusX)
+SaveAs(m, "m_+X_PBC.ovf")
+for iz := 0; iz < 2*N-1; iz++ {
+    for iy := 0; iy < 2*N-1; iy++ {
+        totalCells := 0
+        for ix := 0; ix < 2*N-1; ix++ {
+            totalCells += geom.GetCell(ix, iy, iz)
+        }
+        v := Index2Coord(0, iy, iz)
+        r := Pow(2*Abs(v.Y())/d, 3) + Pow(Abs(2*(Mod(Abs(v.Z()) + N*c/2, N*c) - N*c/2))/d, 3)
+        if ((r <= 1) && (r >= (1 - Pow((N-1)*c/d, 3))) && iy > N/2 && iy < 3*N/2) { // Superball shape in YZ-plane
+            Expect("+X_PBC: One cell thick along X", totalCells, 2, 0)
+        } else {
+            Expect("+X_PBC: No cells along X", totalCells, 0, 0)
+        }
+    }
+}
+
+SetGeom(PlusX.Add(PlusY))
+SaveAs(m, "m_+X+Y.ovf")
+for iz := 0; iz < 2*N-1; iz++ {
+    for ix := 0; ix < 2*N-1; ix++ {
+        totalCells := 0
+        for iy := 0; iy < 2*N-1; iy++ {
+            totalCells += geom.GetCell(ix, iy, iz)
+        }
+        v := Index2Coord(ix, 0, iz)
+        r := Pow(Abs(2*(Mod(Abs(v.X()) + N*c/2, N*c) - N*c/2))/d, 3) + Pow(Abs(2*(Mod(Abs(v.Z()) + N*c/2, N*c) - N*c/2))/d, 3)
+        if (r <= 1 - Pow((N+1)*c/d, 3)) { // Superball shape in XZ-plane
+            Expect("+X+Y_PBC: One cell thick along X", totalCells, 1, 0)
+        } else if (r <= 1) {
+            Expect("+X+Y_PBC: Area with possibly multiple cells along X", sign(totalCells), 1, 0)
+        } else {
+            Expect("+X+Y_PBC: No cells along X", totalCells, 0, 0)
+        }
+    }
+}
+
+SetGeom(PlusX.Add(PlusY.Intersect(MinZ)))
+SaveAs(m, "m_+X+Y-Z.ovf")