-
Notifications
You must be signed in to change notification settings - Fork 149
/
generators.go
323 lines (296 loc) · 8.36 KB
/
generators.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
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
// Copyright 2016 The Gosl Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package msh
import (
"math"
"github.com/cpmech/gosl/chk"
)
// GenQuadRegionHL generates 2D region made of quads (high-level version of GenQuadRegion)
// NOTE: see GenQuadRegion for more details
func GenQuadRegionHL(ctype, ndivR, ndivS int, xmin, xmax, ymin, ymax float64) (o *Mesh) {
return GenQuadRegion(ctype, ndivR, ndivS, false, func(i, j, nr, ns int) (x, y float64) {
dx := (xmax - xmin) / float64(nr-1)
dy := (ymax - ymin) / float64(ns-1)
x = xmin + float64(i)*dx
y = ymin + float64(j)*dy
return
})
}
// GenRing2d generates mesh of quads representing a 2D ring
// ctype -- one of Type{Qua4,Qua8,Qua9,Qua12,Qua16,Qua17}
// ndivR -- number of divisions along radius
// ndivA -- number of divisions along alpha
// r -- minimum radius
// R -- maximum radius
// alpha -- maximum alpha
// NOTE: a circular region is created if maxA=2⋅π
func GenRing2d(ctype int, ndivR, ndivA int, r, R, alpha float64) (o *Mesh) {
circle := math.Abs(alpha-2.0*math.Pi) < 1e-14
return GenQuadRegion(ctype, ndivR, ndivA, circle, func(i, j, nr, ns int) (x, y float64) {
dr := (R - r) / float64(nr-1)
da := alpha / float64(ns-1)
a := float64(j) * da
l := r + float64(i)*dr
x = l * math.Cos(a)
y = l * math.Sin(a)
return
})
}
// GenQuadRegion generates 2D region made of quads
//
// ctype -- one of Type{Qua4,Qua8,Qua9,Qua12,Qua16,Qua17}
// ndivR -- number of divisions (cells) along r (e.g. x)
// ndivS -- number of divisions (cells) along s (e.g. y)
// circle -- connect last row (s=ndivS) with the previous one (s=0)
//
// f(i,j,nr,ns) -- is a function that computes the (x,y) coordinates of grid nodes
// were nr=ndivR+1 and nr=ndivS+1
//
// example (to generate a rectangle):
//
// f := func(i, j, nr, ns int) (x, y float64) {
// dx := (xmax - xmin) / float64(nr-1)
// dy := (ymax - ymin) / float64(ns-1)
// x = xmin + float64(i)*dx
// y = ymin + float64(j)*dy
// return
// }
//
// The boundaries are tagged as below
//
// 34 3 23 30
// @-----@------@ +-----------+
// | | | |
// | | | |
// 4 @ vertices @ 2 40 | edges | 20
// | | | |
// | | | |
// @-----@------@ +-----------+
// 41 1 12 10
//
func GenQuadRegion(ctype, ndivR, ndivS int, circle bool, f func(i, j, nr, ns int) (x, y float64)) (o *Mesh) {
// check
if ndivR < 1 {
chk.Panic("number of divisions along u must be greater than zero\n")
}
if ndivS < 1 {
chk.Panic("number of divisions along v must be greater than zero\n")
}
// type of cell
typekey := TypeIndexToKey[ctype]
// data about shape
nvCell := NumVerts[ctype] // number of vertices of cell
nEdges := len(EdgeLocalVertsD[ctype]) // number of edges of cell
nvEdge := len(EdgeLocalVertsD[ctype][0]) // number of vertices along edge of cell
nvEtot := nvEdge*nEdges - nEdges // total number of vertices along all edges
nvCen := nvCell - nvEtot // number of central vertices of cell
nlCen := nvEdge - 2 // number of lines along central part of cell
nvlCen := make([]int, nlCen) // number of vertices along central line of cell
for i := 0; i < nlCen; i++ {
nvlCen[i] = 2
}
// compute total number of vertices
nvertR := make([]int, 1+nlCen) // number of vertices along r-lines of mesh
nvertR[0] = nvEdge*ndivR - (ndivR - 1) // number of vertices along r-line
for i := 0; i < nlCen; i++ {
nvertR[1+i] = 2*ndivR - (ndivR - 1) // number of vertices along central lines
}
// extra central vertices
if nvCen > 0 {
switch ctype {
case TypeQua9:
nvlCen[0]++ // per cell
nvertR[1] += 1 * ndivR // total along cen line
case TypeQua16:
nvlCen[0] += 2 // per cell
nvlCen[1] += 2 // per cell
nvertR[1] += 2 * ndivR // total along cen line
nvertR[2] += 2 * ndivR // total along cen line
case TypeQua17:
nvlCen[1]++ // per cell
nvertR[2] += 1 * ndivR // total along cen line
default:
chk.Panic("cannot handle central vertices of %q cells\n", typekey)
}
}
// total number of vertices
nvTot := 0
for _, nr := range nvertR {
nvTot += nr * ndivS
}
if !circle {
nvTot += nvertR[0] // verts along top "main" line
}
// total number of lines along s
ns := (1+nlCen)*ndivS + 1 // +1 => top line (needed even with circle)
// allocate vertices
o = new(Mesh)
o.Verts = make([]*Vertex, nvTot)
// generate vertices
j, iv := 0, 0
for j < ns-1 {
for _, nr := range nvertR {
for i := 0; i < nr; i++ {
x, y := f(i, j, nr, ns)
o.Verts[iv] = &Vertex{ID: iv, Tag: qvtag(i, j, nr, ns), X: []float64{x, y}}
iv++
}
j++
}
}
if !circle {
nr := nvertR[0]
for i := 0; i < nr; i++ {
x, y := f(i, ns-1, nr, ns)
o.Verts[iv] = &Vertex{ID: iv, Tag: qvtag(i, ns-1, nr, ns), X: []float64{x, y}}
iv++
}
}
// allocate cells
ncTot := ndivR * ndivS
o.Cells = make([]*Cell, ncTot)
// constants
sumNvr := 0
for _, nr := range nvertR {
sumNvr += nr
}
// generate cells
ic := 0
for j := 0; j < ndivS; j++ {
for i := 0; i < ndivR; i++ {
// compute pivot points
incR := i * (nvEdge - 1) // increment along r
remR := nvertR[0] - incR // remainder along r +1
a := incR + sumNvr*j // lower pivot
b := incR + sumNvr*(j+1) // upper pivot
up := func(idxCenLine int) (idx int) {
idx = a + remR
for n := idxCenLine - 1; n >= 0; n-- {
idx += nvertR[1+n]
}
idx += (nvlCen[idxCenLine] - 1) * i
return
}
// correct upper pivot
if circle && j == ndivS-1 {
b = incR
}
// set connectivity
var v []int
switch ctype {
case TypeQua4:
v = []int{a, a + 1, b + 1, b}
case TypeQua8:
c := up(0)
v = []int{a, a + 2, b + 2, b, a + 1, c + 1, b + 1, c}
case TypeQua9:
c := up(0)
v = []int{a, a + 2, b + 2, b, a + 1, c + 2, b + 1, c, c + 1}
case TypeQua12:
c := up(0)
d := up(1)
v = []int{a, a + 3, b + 3, b, a + 1, c + 1, b + 2, d, a + 2, d + 1, b + 1, c}
case TypeQua16:
c := up(0)
d := up(1)
v = []int{a, a + 3, b + 3, b, a + 1, c + 3, b + 2, d, a + 2, d + 3, b + 1, c, c + 1, c + 2, d + 2, d + 1}
case TypeQua17:
c := up(0)
d := up(1)
e := up(2)
v = []int{a, a + 4, b + 4, b, a + 1, c + 1, b + 3, e, a + 2, d + 2, b + 2, d, a + 3, e + 1, b + 1, c, d + 1}
default:
chk.Panic("cannot handle cell type = %q at this time\n", typekey)
}
// set cell
o.Cells[ic] = &Cell{ID: ic, Tag: -1, TypeKey: typekey, V: v, EdgeTags: qetag(i, j, ndivR, ndivS)}
ic++
}
}
// results
o.CheckAndCalcDerivedVars()
return
}
// auxiliary //////////////////////////////////////////////////////////////////////////////////////
// qvtag returns vertex tags for quad-meshes
//
// example:
//
// 34 3 3 3 23
// @---@---@---@---@
// | |
// 4 @ @ 2
// | |
// 4 @ @ 2
// | |
// 4 @ @ 2
// | |
// @---@---@---@---@
// 41 1 1 1 12
//
// nx and ny are number of lines (rows/cols)
func qvtag(i, j, nx, ny int) int {
if i == 0 {
if j == 0 {
return 41
}
if j == ny-1 {
return 34
}
return 4
}
if j == 0 {
if i == nx-1 {
return 12
}
return 1
}
if i == nx-1 {
if j == ny-1 {
return 23
}
return 2
}
if j == ny-1 {
return 3
}
return 0
}
// qetag returns edge tags for quad-meshes
//
// example:
//
// 30
// +------------+
// | |
// | |
// 40 | | 20
// | |
// | |
// +------------+
// 10
//
// nx and ny are number of cells along each direction
func qetag(i, j, nx, ny int) []int {
if nx*ny == 1 {
return []int{10, 20, 30, 40}
}
if i > 0 && i < nx-1 && j > 0 && j < ny-1 {
return nil
}
tags := []int{0, 0, 0, 0}
if j == 0 {
tags[0] = 10
}
if i == nx-1 {
tags[1] = 20
}
if j == ny-1 {
tags[2] = 30
}
if i == 0 {
tags[3] = 40
}
return tags
}