-
Notifications
You must be signed in to change notification settings - Fork 0
/
ext_make3dgrains.go
159 lines (133 loc) · 4.6 KB
/
ext_make3dgrains.go
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
// 3D Voronoi tessellation. Contributed by Peyton Murray.
package engine
import (
"math"
"math/rand"
)
func init() {
DeclFunc("ext_make3dgrains", Voronoi3d, "3D Voronoi tesselation over shape (grain size, shellthickness, gbthickness, starting region number, num regions, shape, seed)")
}
func Voronoi3d(grainsize float64, shellthickness float64, gbthickness float64, startRegion int, numRegions int, inputShape Shape, seed int) {
Refer("Lel2014")
SetBusy(true)
defer SetBusy(false)
t := newTesselation3d(grainsize, numRegions, int64(seed), startRegion, inputShape)
regions.hist = append(regions.hist, t.RegionOf)
regions.render(t.RegionOf)
}
type tesselation3d struct {
grainsize float64
maxRegion int
rnd *rand.Rand
startRegion int
shape Shape //Shape of the tesselated region
centers []center3d //List of Voronoi centers
}
// Stores location of each Voronoi center
type center3d struct {
x, y, z float64 // center position (m)
region byte // region for all cells near center
}
// Stores location of each cell
type cellLocs struct{ x, y, z float64 }
// nRegion exclusive
func newTesselation3d(grainsize float64, nRegion int, seed int64, startRegion int, inputShape Shape) *tesselation3d {
t := tesselation3d{grainsize,
nRegion,
rand.New(rand.NewSource(seed)),
startRegion,
inputShape,
make([]center3d, 0)}
t.makeRandomCenters()
return &t
}
// Permutes the slice of cell locations. I don't understand why this needs to be done if we're choosing
// random (Intn()) cells out of the slice of cell locations, but hey, it seems to do the trick.
func shuffleCells(src []cellLocs) []cellLocs {
dest := make([]cellLocs, len(src))
perm := rand.Perm(len(src))
for i, v := range perm {
dest[v] = src[i]
}
return dest
}
func (t *tesselation3d) makeRandomCenters() {
//Make a list of all the cells in the shape.
cells := t.tabulateCells()
cells = shuffleCells(cells)
//Choose number of grains to make. Assume volume of grain is given by (4/3)*pi*r^3
shapeVolume := cellVolume() * float64(len(cells))
grainVolume := (float64(1) / 6) * math.Pi * t.grainsize * t.grainsize * t.grainsize
nAvgGrains := shapeVolume / grainVolume
nGrains := t.truncNorm(nAvgGrains)
//TODO: same cell can be chosen twice by random chance
t.centers = make([]center3d, nGrains)
for p := 0; p < nGrains; p++ {
rndCell := cells[t.rnd.Intn(nGrains)]
t.centers[p].x = rndCell.x
t.centers[p].y = rndCell.y
t.centers[p].z = rndCell.z
randRegion := t.startRegion + t.rnd.Intn(t.maxRegion)
t.centers[p].region = byte(randRegion)
}
return
}
// Creates a slice of all cells which fall in the shape specified in the constructor.
func (t *tesselation3d) tabulateCells() []cellLocs {
//Initialze array of cells
cells := make([]cellLocs, 0)
//Get the mesh size
meshSize := MeshSize()
//Iterate across all cells in the mesh, and append those that are inside the shape
for ix := 0; ix < meshSize[0]; ix++ {
for iy := 0; iy < meshSize[1]; iy++ {
for iz := 0; iz < meshSize[2]; iz++ {
cell := Index2Coord(ix, iy, iz)
x := cell.X()
y := cell.Y()
z := cell.Z()
if t.shape(x, y, z) {
cells = append(cells, cellLocs{x, y, z})
}
}
}
}
print("Number of cells in region: ", len(cells), "\n")
print("Number of cells in universe: ", meshSize[0]*meshSize[1]*meshSize[2], "\n")
return cells
}
// Find the nearest Voronoi center to the point (x, y, z). Only points inside the given shape will be
// assigned a region.
func (t *tesselation3d) RegionOf(x, y, z float64) int {
if t.shape(x, y, z) {
nearest := center3d{x, y, z, 0}
mindist1 := math.Inf(1)
mindist2 := math.Inf(1)
for _, c := range t.centers {
dist := sqrt(sqr(x-c.x) + sqr(y-c.y) + sqr(z-c.z))
if dist < mindist1 {
nearest = c
mindist2 = mindist1
mindist1 = dist
} else if dist < mindist2 {
mindist2 = dist
}
}
sindex := math.Floor(math.Cos(math.Floor((mindist2-mindist1)/(2*shellthickness+gbthickness))))
bindex := math.Floor(math.Cos(math.Floor((mindist2-mindist1)/gbthickness)))
return int(nearest.region+(sindex+bindex)*numRegions)
} else {
return -1 //When the regions are rendered, any region < 0 will not be rastered.
}
}
// Generate normally distributed numbers; mean = lambda, variance = lambda. If generated number < 0, return 1.
// Equivalent to Poisson distribution (with mean = lambda) for large lambda (which is usually true, since the volume
// of a grain is usually much less than the simulation volume.
func (t *tesselation3d) truncNorm(lambda float64) int {
ret := lambda + math.Sqrt(lambda)*t.rnd.NormFloat64()
if ret <= 0 {
return 1
} else {
return int(ret + 0.5)
}
}