Line data Source code
1 : module Meshes
2 :
3 : import PMesh: Mesh, Types, VHnd, nv, nf, vattr, addvertex!, addvertices!,
4 : addface!, insertedge!, halfedge, istrianglemesh, createmesh, vertices
5 :
6 : const TEMPLATE = Types.MESH3N
7 : const V3 = Types.V3
8 :
9 : const VQUAD = [V3(0, 0, 0), V3(1, 0, 0), V3(0, 1, 0), V3(1, 1, 0)]
10 :
11 : """
12 : mesh = triangle([; template::Mesh])
13 :
14 : Create planar 3d mesh with one triangle.
15 : """
16 14 : function triangle(; template::Mesh = TEMPLATE)
17 7 : mesh = similar(template)
18 7 : a, b, c = VHnd.(addvertices!(mesh, VQUAD[1:3]))
19 7 : addface!(mesh, a, b, c)
20 :
21 7 : @assert nv(mesh) == 3 && nf(mesh) == 1
22 :
23 7 : mesh
24 : end
25 :
26 : """
27 : mesh = quad([; template::Mesh])
28 :
29 : Create planar 3d mesh with one quad.
30 : """
31 6 : function quad(; template::Mesh = TEMPLATE)
32 3 : mesh = similar(template)
33 3 : a, b, c, d = VHnd.(addvertices!(mesh, VQUAD))
34 3 : addface!(mesh, a, b, c, d)
35 :
36 3 : @assert nv(mesh) == 4 && nf(mesh) == 1
37 :
38 3 : mesh
39 : end
40 :
41 : """
42 : mesh = wing([; template::Mesh])
43 :
44 : Create planar 3d mesh with two triangles, which triangulate a quad.
45 : """
46 2 : function wing(; template::Mesh = TEMPLATE)
47 1 : mesh = similar(template)
48 1 : a, b, c, d = VHnd.(addvertices!(mesh, VQUAD))
49 1 : addface!(mesh, a, b, c)
50 1 : addface!(mesh, c, b, d)
51 :
52 1 : @assert nv(mesh) == 4 && nf(mesh) == 2
53 :
54 1 : mesh
55 : end
56 :
57 : """
58 : mesh = nring([;n = 6, r=A(t), α=B(t), template::Mesh])
59 :
60 : Create 3d mesh of a 1-ring with `n` neighbors around
61 : the center vertex `VHnd(1)` at `(0, 0, 0)`.
62 :
63 : The functions `r(t)` and `α(t)` determine the vertex positions
64 : for `t = [0, 1/n, ..., (n-1)/m]`. The default positions boundary
65 : vertices regularly on the unit circle in the `x-y-plane`.
66 :
67 : See also [`ngon`](@ref)
68 : """
69 32 : function nring(;n::Int=6, r=t -> 1, α=t ->2π*t, template::Mesh = TEMPLATE)
70 3 : @assert n > 1
71 3 : mesh = similar(template)
72 3 : x = vattr(mesh, Val(:x))
73 3 : v = addvertex!(mesh)
74 3 : x[v] = V3(0, 0, 0)
75 :
76 3 : vs = [addvertex!(mesh) for _ in 1:n]
77 :
78 3 : for i in 1:n
79 26 : t = (i - 1) / n
80 26 : rt = r(t)
81 26 : αt = α(t)
82 26 : x[vs[i]] = V3(rt * cos(αt), rt * sin(αt), 0)
83 :
84 26 : j = (i % n) + 1
85 26 : addface!(mesh, v, vs[i], vs[j])
86 49 : end
87 :
88 3 : @assert nv(mesh) == n+1 && nf(mesh) == n
89 :
90 3 : mesh
91 : end
92 :
93 : """
94 : mesh = ngon([;n = 6, r=A(t), α=B(t), template::Mesh])
95 :
96 : Create 3d mesh of a `n`-gon.
97 :
98 : The functions `r(t)` and `α`(t)` determine the vertex positions
99 : for `t = [0, 1/n, ..., (n-1)/m]`. The default positions boundary
100 : vertices regularly on the unit circle in the `x-y-plane`.
101 :
102 : See also [`nring`](@ref)
103 : """
104 24 : function ngon(;n::Int=6, r=t -> 1, α=t ->2π*t, template::Mesh = TEMPLATE)
105 2 : @assert n > 1
106 2 : mesh = similar(template)
107 2 : x = vattr(mesh, Val(:x))
108 :
109 2 : vs = [addvertex!(mesh) for _ in 1:n]
110 :
111 2 : for i in 1:n
112 20 : t = (i - 1) / n
113 20 : rt = r(t)
114 20 : αt = α(t)
115 20 : x[vs[i]] = V3(rt * cos(αt), rt * sin(αt), 0)
116 38 : end
117 :
118 2 : addface!(mesh, vs)
119 :
120 2 : @assert nv(mesh) == n && nf(mesh) == 1
121 :
122 2 : mesh
123 : end
124 :
125 : """
126 : mesh = cube([; triangulate=false, template::Mesh])
127 :
128 : Create standard cube with each quad optionally split into two `triangles`.
129 : """
130 6 : function cube(;triangles=false, template::Mesh = TEMPLATE)
131 3 : mesh = similar(template)
132 3 : x = vattr(mesh, Val(:x))
133 :
134 3 : vs = [addvertex!(mesh) for _ in 1:8]
135 6 : x[vs[1:4]] = VQUAD[[1, 2, 4, 3]]
136 18 : x[vs[5:8]] = map(x -> V3(x[1], x[2], 1), VQUAD[[1, 2, 4, 3]])
137 :
138 3 : addface!(mesh, vs[1], vs[2], vs[3], vs[4])
139 3 : addface!(mesh, vs[5], vs[8], vs[7], vs[6])
140 :
141 3 : addface!(mesh, vs[1], vs[4], vs[8], vs[5])
142 3 : addface!(mesh, vs[2], vs[6], vs[7], vs[3])
143 :
144 3 : addface!(mesh, vs[1], vs[5], vs[6], vs[2])
145 3 : addface!(mesh, vs[4], vs[3], vs[7], vs[8])
146 :
147 3 : @assert nv(mesh) == 8 && nf(mesh) == 6
148 :
149 3 : if triangles
150 6 : SPLIT = [ ((VHnd(1), VHnd(2)), (VHnd(3), VHnd(4))),
151 : ((VHnd(6), VHnd(5)), (VHnd(8), VHnd(7))),
152 : ((VHnd(1), VHnd(4)), (VHnd(8), VHnd(5))),
153 : ((VHnd(3), VHnd(2)), (VHnd(6), VHnd(7))),
154 : ((VHnd(1), VHnd(5)), (VHnd(6), VHnd(2))),
155 : ((VHnd(8), VHnd(4)), (VHnd(3), VHnd(7)))]
156 :
157 :
158 1 : for (a, b) in SPLIT
159 6 : insertedge!(mesh, halfedge(mesh, a), halfedge(mesh, b))
160 6 : end
161 :
162 1 : @assert nv(mesh) == 8 && nf(mesh) == 12 && istrianglemesh(mesh)
163 : end
164 :
165 3 : mesh
166 : end
167 :
168 : """
169 : mesh = tetrahedron([; template::Mesh])
170 :
171 : Create standard tetrahedron.
172 : """
173 4 : function tetrahedron(; template::Mesh = TEMPLATE)
174 2 : mesh = similar(template)
175 2 : x = vattr(mesh, Val(:x))
176 2 : a, b, c, d = VHnd.(addvertices!(mesh, VQUAD))
177 2 : x[d] = V3(0, 0, 1)
178 :
179 2 : addface!(mesh, a, b, c)
180 2 : addface!(mesh, b, a, d)
181 2 : addface!(mesh, c, b, d)
182 2 : addface!(mesh, a, c, d)
183 :
184 2 : @assert nv(mesh) == 4 && nf(mesh) == 4
185 :
186 2 : mesh
187 : end
188 :
189 : #------------------------------------------------------------------------------
190 :
191 : #
192 : # code from lecture mp23 (solution to assignment)
193 : #
194 :
195 3 : function _rowsperiodic!(idx)
196 3 : idx[end, :] .= idx[1, :]
197 3 : idx
198 : end
199 :
200 14 : function _colsperiodic!(idx)
201 14 : idx[:, end] .= idx[:, 1]
202 14 : idx
203 : end
204 :
205 8 : function _rowpoles!(idx)
206 128 : idx[1, :] .= idx[1]
207 128 : idx[end, :] .= idx[end]
208 8 : idx
209 : end
210 :
211 10176 : _validtriangle(i, j, k) = (i != j && j != k && k != i)
212 :
213 : """
214 : x, t = _samplemesh(f, u, v; boundary=identity)
215 :
216 : Sample (x, y, z) = f(u, v)` given ranges `u` and `v`, apply
217 : tessellation and return mesh `x, t` in shared vertex representation.
218 :
219 : The function `boundary` reassigns vertex indices: `_rowsperiodic` and
220 : `_colsperiodic!` create periodic "boundaries", `_rowpoles!` maps `v`
221 : boundaries to a single vertex. The `boudary` function takes an index
222 : grid as input and returns the same grid to enable concatenation, e.g.,
223 : `_rowsperiodic! ∘ _colsperiodic!` for a torus.
224 :
225 : The method "filters: unused vertex indices and invalid triangles.
226 : """
227 32 : function _samplemesh(f, u, v; boundary=identity)
228 16 : n, m = length(u), length(v)
229 16 : T = eltype(u)
230 :
231 16 : idxgrid = collect(reshape(1:(m * n), (m, n)))
232 :
233 : # apply boundary configuration
234 16 : boundary(idxgrid)
235 :
236 : # filter unused indices: "re-index"
237 5476 : idxused = fill(false, length(idxgrid))
238 5476 : idxused[vec(idxgrid)] .= true
239 16 : nused = count(idxused)
240 5476 : idxmap = zeros(Int, length(idxgrid))
241 32 : idxmap[idxused] .= 1:nused
242 5492 : idxgrid = map(idx -> idxmap[idx], idxgrid)
243 :
244 : # sample surface f(u, v)
245 15126 : x = zeros(T, 3, nused)
246 :
247 16 : ci = 1
248 16 : for (k, idx) in enumerate(CartesianIndices(idxgrid))
249 5476 : idxused[k] || continue
250 5042 : i, j = Tuple(idx)
251 :
252 5188 : x[:, ci] .= f(u[j], v[i])
253 5042 : ci += 1
254 10936 : end
255 :
256 : # tessellate regular mesh
257 30528 : t = zeros(Int, 3, 2 * (m - 1) * (n - 1))
258 :
259 16 : ci = 0
260 :
261 10192 : function _add(a, b, c)
262 10416 : _validtriangle(a, b, c) || return
263 :
264 9936 : ci += 1
265 9936 : t[:, ci] .= (a, b, c)
266 : end
267 :
268 16 : for j in 1:(n-1), i in 1:(m-1)
269 5088 : _add(idxgrid[i, j], idxgrid[i+1, j], idxgrid[i, j+1])
270 5088 : _add(idxgrid[i+1, j], idxgrid[i+1, j+1], idxgrid[i, j+1])
271 5256 : end
272 :
273 16 : x, copy(t[:, 1:ci])
274 : end
275 :
276 4324 : function _sphere(r::T, u::T, v::T) where {T<:Real}
277 4324 : (r * sin(v) * cos(u),
278 : r * sin(v) * sin(u),
279 : r * cos(v))
280 : end
281 :
282 : """
283 : mesh = sphere([; r=1, m=8 ,n=m, triangles=true, template::Mesh])
284 :
285 : Generate sphere mesh with radius `r` tessellated on a `m × n`
286 : grid of faces (with two poles).
287 :
288 : See also [`torus`](@ref)
289 : """
290 12 : function sphere(; r=1.0, m::Int=8, n::Int=m, triangles::Bool=true,
291 : template::Mesh = TEMPLATE)
292 6 : u = range(0, 2π, length=n + 1)
293 6 : v = range(0, π, length=m + 1)
294 :
295 4330 : f = (u, v) -> _sphere(r, u, v) # bind radius
296 6 : _samplemesh(f, u, v; boundary = _rowpoles! ∘ _colsperiodic!) |>
297 : _tessellate(triangles; template)
298 : end
299 :
300 192 : function _torus(R::T, r::T, u::T, v::T) where {T<:Real}
301 192 : (( R + r * cos(v)) * cos(u),
302 : ( R + r * cos(v)) * sin(u),
303 : -r * sin(v))
304 : end
305 :
306 : """
307 : mesh = torus([; ro=1, ri=0.25, m=8, n=8, triangles=true, template::Mesh])
308 :
309 : Generate torus mesh with outer radius `ro and inner ("tube")
310 : radius `ri` tessellated on a `m × n` grid of faces.
311 :
312 : See also [`sphere`](@ref)
313 : """
314 6 : function torus(; ro=1.0, ri=0.25, m::Int=8, n::Int=m, triangles::Bool=true,
315 : template::Mesh = TEMPLATE)
316 3 : @assert ro >= ri
317 :
318 3 : u = range(0, 2π, length=n + 1)
319 3 : v = range(0, 2π, length=m + 1)
320 :
321 195 : f = (u, v) -> _torus(ro, ri, u, v) # bind radii
322 3 : _samplemesh(f, u, v;
323 : boundary=_rowsperiodic! ∘ _colsperiodic!) |>
324 : _tessellate(triangles; template)
325 : end
326 :
327 : """
328 : mesh = tessellate(x, t[; template::Mesh])
329 :
330 : Create mesh from shared vertex table using `TEMPLATE`.
331 :
332 : See also [`createmesh`](@ref)
333 : """
334 32 : function tessellate(x, t; template = TEMPLATE)
335 16 : mesh = createmesh(t, template)
336 :
337 16 : mx = vattr(mesh, Val(:x))
338 16 : Ns = 1:length(mx[VHnd(1)])
339 :
340 32 : for (j, v) in enumerate(vertices(mesh))
341 5042 : mx[v] = V3(x[:, j])[Ns]
342 10068 : end
343 :
344 16 : mesh
345 : end
346 :
347 38 : function _tessellate(triangles; template)
348 16 : if triangles
349 20 : xt -> tessellate(xt...; template)
350 : else
351 12 : function f(xt)
352 6 : x, t = xt
353 6 : m = size(t, 2)
354 3008 : qs = zeros(Int, 4, m)
355 :
356 6 : j = 1
357 6 : k = 1
358 403 : while j < m
359 1191 : ta = t[:, j]
360 1191 : tb = t[:, j+1]
361 :
362 : # guess triangulated quad from _samplemesh
363 397 : if ta[2] == tb[1] && ta[3] == tb[3]
364 352 : qs[:, k] .= (t[1, j], t[2, j], t[2, j+1], t[3, j])
365 352 : j += 2
366 : else
367 45 : qs[1:3, k] .= t[:, j]
368 45 : j += 1
369 : end
370 :
371 397 : k += 1
372 397 : end
373 :
374 6 : if j == m
375 3 : qs[1:3, k] = t[:, j]
376 3 : k += 1
377 : end
378 :
379 6 : tessellate(x, qs[:, 1:k-1]; template)
380 : end
381 :
382 : end
383 : end
384 :
385 360 : function _cylinder(r::T, u::T, v::T) where {T<:Real}
386 360 : (r * cos(u),
387 : r * sin(u),
388 : v)
389 : end
390 :
391 148 : function _cylinderclosed(r::T, h::T, u::T, v::T) where {T<:Real}
392 148 : (v < 0) && return (T(0), T(0), T(0))
393 146 : (v > h) && return (T(0), T(0), h)
394 :
395 144 : _cylinder(r, u, v)
396 : end
397 :
398 : """
399 : mesh = cylinder([; rh=1, h=1, m=8, n=8,
400 : triangles=true, closed=false, template::Mesh])
401 :
402 : Generate cylinder mesh with radius `r and height `h`
403 : tessellated on a `m × n` grid of faces.
404 :
405 : See also [`sphere`](@ref)
406 : """
407 10 : function cylinder(; r=1.0, h=1.0, m::Int=8, n::Int=m,
408 : triangles::Bool=true, closed=false, template::Mesh = TEMPLATE)
409 5 : closed && return _closedcylinder(; r, h, m, n, triangles, template)
410 :
411 3 : _opencylinder(; r, h, m, n, triangles, template)
412 : end
413 :
414 4 : function _closedcylinder(; r, h, m, n, triangles::Bool=true, template)
415 2 : u = range(0, 2π, length=n + 1)
416 2 : v = collect(range(0, h, length=m + 1))
417 2 : v = vcat([-h], v, [h+1])
418 :
419 296 : f = (u, v) -> _cylinderclosed(r, h, u, v) # bind radius, h
420 2 : _samplemesh(f, u, v; boundary = _rowpoles! ∘ _colsperiodic!) |>
421 : _tessellate(triangles; template)
422 : end
423 :
424 6 : function _opencylinder(; r, h, m, n, triangles::Bool=true, template)
425 3 : u = range(0, 2π, length=n + 1)
426 3 : v = range(0, h, length=m + 1)
427 :
428 219 : f = (u, v) -> _cylinder(r, u, v) # bind radius
429 3 : _samplemesh(f, u, v; boundary = _colsperiodic!) |>
430 : _tessellate(triangles; template)
431 : end
432 :
433 : """
434 : mesh = rectangle([; w=1.0, h=1.0, triangles=true, z=(x, y) -> 0.0,
435 : template::Mesh])
436 :
437 : Generate a rectangular patch with width `w` and height `h` tessellated
438 : on a `m × n` grid of faces. Optionally lift third coordinate to
439 : `z(x, y)`.
440 :
441 : See also [`sphere`](@ref)
442 : """
443 4 : function rectangle(; w=1.0, h=1.0, m::Int=8, n::Int=m, triangles=true,
444 81 : z = (u, v) -> zero(u), template::Mesh = TEMPLATE)
445 2 : u = range(0, w, length=n + 1)
446 2 : v = range(0, h, length=m + 1)
447 :
448 164 : f = (u, v) -> (u, v, z(u, v))
449 2 : _samplemesh(f, u, v) |> _tessellate(triangles; template)
450 : end
451 :
452 : # TODO: right w/ hole (2 boundaries)
453 :
454 : end
|
